{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 統計的推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.320204Z",
     "start_time": "2018-08-14T18:04:46.813094Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "%precision 3\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.327204Z",
     "start_time": "2018-08-14T18:04:47.321803Z"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('../data/ch4_scores400.csv')\n",
    "scores = np.array(df['点数'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.332144Z",
     "start_time": "2018-08-14T18:04:47.328296Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(69.530, 206.669)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_mean = np.mean(scores)\n",
    "p_var = np.var(scores)\n",
    "\n",
    "p_mean, p_var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.534295Z",
     "start_time": "2018-08-14T18:04:47.333230Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAFoCAYAAABpKJ6VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl0lPed5/v3tzZJgJAMkowkhMRiE4PNKhbbGLCxcRhvjeks7sR9k44HZ266z73dJ33b03NnznR62smZzNy56c6cdnD3jdMdO3HsjsdJZGwDBmyw2TeD2Y3EJjBISICQVNvv/lEFjaEEpa2ekvR5ncNBVfzqqc+jH1V8eOrR7zHnHCIiIiKSGT6vA4iIiIgMJCpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhkU8DrAjRQVFbmqqqpefY5wOEwoFOrV55DO07xkH81JdtK8ZB/NSXbKxLxs3br1rHOu+Gbjsrp8VVVVsWXLll59jtraWnq74EnnaV6yj+YkO2leso/mJDtlYl7MrC6dcfrYUURERCSDVL5EREREMkjlS0RERCSDVL5EREREMkjlS0RERCSDVL5EREREMkjlS0RERCSDVL5EREREMkjlS0RERCSDblq+LOFFM9tmZt/rYEyOmb1hZtvN7Nlr/uzHZvbL5NfDzGyVme0ws8d6ZhdERERE+o50jnzNBnKBamCpmZWlGPMksBuYDzxvZjkAZvYQMPSqcc8CrwBfB37Y9dgiIiIifVM65Ws6sBqYBKwDpt5gzGxgAzDOzAqB7wB/nWLcBKDRzAZ3PbqIiIhI35NO+SoAzgIlwJrk7Y7GOGBv8vYPgeeASIpx54CTHWxLREREpN8KpDGmGcgBPiFx5OpIB2OKgJ3AEqAFmAPcRuIjyzFm9nRy3FgSR9C+m7z9OWa2FFgKUF5eTm1tbad2qLMaGhp6dfvSNZqX7KM5yU7ZPi/zX9hz3X1rvj3RgySZk+1zMlBl07ykU762Akucc6+Z2VzgZx2MmUDiyNgkYJ9z7g4AM6sCfuCc+2czGwmUAweAfOdcy7Ubcs4tA5YBVFdXu6qqqk7uUudl4jmk8zQv2Udzkp2ye16uL1/ZnbdnDIR97IuyZV7S+dhxA1BoZhuBzc65kynGvAEsSI59yTnX3sG2lgF/SqKkfb/zcUVERET6tpse+XLOOeCZy7fNbB6wChjrnKtLjmkDFnfw+Frgq8mvG0iUNBEREZEBqSuLrG4BppA4YV5EREREOiGdc74+J3me1u5eyCIiIiLS7+nyQiIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZpPIlIiIikkEqXyIiIiIZdNPyZQkvmtk2M/teB2NyzOwNM9tuZs8m75tsZh+a2Toz+1lyO4+Y2YHkfevMbEhP75CIiIhINkvnyNdsIBeoBpaaWVmKMU8Cu4H5wPNmlgP4gSedc3OAUuBWYAjwvHNuTvLXxR7YBxEREZE+I53yNR1YDUwC1gFTbzBmNrABGOec2wbMM7ODwDbn3CkS5euPkkfE/qEndkBERESkL0mnfBUAZ4ESYE3ydkdjHLD38hjn3KvAROALZlYJ1AE/cs7dAwTN7L7u7oCIiIhIXxJIY0wzkAN8QuII15EOxhQBO4ElQLOZTQN2OefCZrYWmJ0sY5dtBMYCH1y9ITNbCiwFKC8vp7a2tlM71FkNDQ29un3pGs1L9tGcZKfuzMv8F/Zcd9+ab0/sTpy09Pb7emf19PdBr5XslE3zkk752goscc69ZmZzgZ91MGYCiSNjk4BDwI+AV0l8HDkVWG9mfwHsds7VkPiI8u+v3ZBzbhmwDKC6utpVVVV1cpc6LxPPIZ2neck+mpPs1PV5ub509PwcZ+I5uqvnM2bfPgpkz7yk87HjBqDQzDYCm51zJ1OMeQNYkBz7knOuHfgvwF+a2Xqg3jm3Cfj/gD8zsy1As3Puox7ZCxEREZE+4qZHvpxzDnjm8m0zmwesAsY65+qSY9qAxdc87jjw0DX3nSFR0kREREQGpK4ssroFmAKkOgImIiIiIjeQzjlfn+OcayGxppeIiIiIdJIuLyQiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQSpfIiIiIhmk8iUiIiKSQQGvA4iISM+qeq7muvtqf/CIJ8/bn3j1fZX+R0e+RERERDJI5UtEREQkg1S+RERERDJI5UtEREQkg1S+RERERDJI5UtEREQkg1S+RERERDJI5UtEREQkg1S+RERERDJI5UtEREQkg1S+RERERDJI5UtEREQkg1S+RERERDJI5UtEREQkg1S+RERERDLopuXLEl40s21m9r0OxuSY2Rtmtt3Mnk3eN9nMPjSzdWb2s+R2rhsnIiIiMpCkc+RrNpALVANLzawsxZgngd3AfOB5M8sB/MCTzrk5QClwawfjRERERAaMdMrXdGA1MAlYB0y9wZjZwAZgnHNuGzDPzA4C25xzp1KN6/YeiIiIiPQh6ZSvAuAsUAKsSd7uaIwD9l4e45x7FZgIfMHMKjsaJyIiIjJQBNIY0wzkAJ+QOHJ1pIMxRcBOYAnQbGbTgF3OubCZrSVxtOu6cdduyMyWAksBysvLqa2t7eQudU5DQ0Ovbl+6RvOSfTQn2SndeUn3vbS333Mz9Rzd1Z3vl14r2Smb5iWd8rUVWOKce83M5gI/62DMBBJHxiYBh4AfAa+S+JhxKrCexJG2a8d9jnNuGbAMoLq62lVVVXVmf7okE88hnad5yT6ak+x0/bzsSWNMZ8al6/rtpZJ9f496/vuVffsokD3zks7HjhuAQjPbCGx2zp1MMeYNYEFy7EvOuXbgvwB/aWbrgXrn3KYOxomIiIgMGDc98uWcc8Azl2+b2TxgFTDWOVeXHNMGLL7mcceBh66577pxIiIiIgNJVxZZ3QJMAVIdARMRERGRG0jnnK/Pcc61kFirS0REREQ6SZcXEhEREckglS8RERGRDFL5EhEREckglS8RERGRDFL5EhEREckglS8RERGRDFL5EhEREckglS8RERGRDFL5EhEREckglS8RERGRDOr05YVERES6xjHYwhw4cIBYLIZzDuccAGZGcXExw4cPx+fTcQHp31S+RESkV/iIM8J3gWJfC0W+Fop9LeRZlF/84uMOHxMMBhkxYgSlpaWUlZVx++23k5eXl8HUIr1P5UtERHrUYGtnvP8MtwfOkmdRnIMml8vxWAFn4oP5h28/SDAYxMyAxFGvWCzG6dOnqa+vp76+nu3bt7Np0yYCgQATJ05k2rRpVFRUXHmMSF+m8iUiIt3mnKPc18wXAp9R4WvGAcfihRyIFnE6nk8E/5Wx9/3PXSm3UfuDR5gyZQoA8XicU6dOsX37dnbt2sXOnTspKSlh2rRpTJ06lVAolIndEukVKl8iItItp06doqamhoU5x2l1AXZFS9kfK6LF5XR5mz6fj7KyMsrKynjooYfYvXs3W7du5e233+bDDz9k4cKFTJgwQUfCpE9S+RIRkS5pb29n9erVbNq0iby8PNaFKzkcG068h3+QPhQKMW3aNKZNm0ZdXR3Lly/n9ddfZ/To0SxatIji4uIefT6R3qYfKRERkU5xzrF7925+/OMfs3HjRqZPn84f//EfczBW3OPF61qVlZUsXbqURYsWUV9fzwsvvMC7775LOBzu1ecV6Uk68iUiImlrb2/nzTffZO/evZSWlvLVr36V8vLyjGbw+XzMnDmTiRMnsnLlSj766CMOHjzIV77yFYqKijKaRaQrdORLRETS0tDQwD/+4z+yb98+HnzwQZ555pmMF6+rDR48mCeeeIKnn36aS5cu8eKLL7J3717P8oikS+VLRERu6sCBA7z44ou0tLTw9NNPc++992bNYqhjxoxh6dKlFBcX86tf/YqVK1cSj8e9jiXSIX3sKCIiHXLO8f7777NmzRpKS0v58pe/TGFhodexrlNQUMA3vvENli9fzvr166mvr2fJkiUMGjTI62gi18mO/7aIiEjWicVi/PrXv2bNmjVMnjyZb37zm1lZvC4LBAI89thjPP7449TV1fHTn/6UCxcueB1L5DoqXyIicp1YLMbrr7/O7t27efDBB3niiScIBoNex0rL1KlTefrppzl//jw//elPaWpq8jqSyOeofImIyOdEo1FeffVV9u3bxxe/+EXuvffePreYaWVlJU8//TStra289NJLNDY2eh1J5AqVLxERuSISifDLX/6SgwcP8uijjzJr1iyvI3XZyJEj+cM//EPC4TA//elPOXPmjNeRRACVLxERSQoQ45VXXuHw4cM8/vjjTJ8+3etI3VZaWso3vvENnHO89NJLnDp1yutIIipfIiICPuIsCB2irq6OxYsXM3XqVK8j9ZiSkhK++c1vEggE+PnPf865c+e8jiQDnMqXiMiA55gTrKXMf4EnnniCSZMmeR2oxw0fPpyvf/3rxGIxXn75ZS5duuR1JBnAVL5ERAa46YETjA00siVSzuTJk72O02uKi4t56qmnaGpq4he/+AWRSMTrSDJA3bR8WcKLZrbNzL7XwZgcM3vDzLab2bPJ+0ab2Soz22Jmf5W87ztm9rGZrTOz5T27KyIi0lnj/Z8xKXiKfdFiPo6O8DpOrxs1ahRLlizh+PHj/Mu//ItWwhdPpHPkazaQC1QDS82sLMWYJ4HdwHzgeTPLAf4D8O+dc9XAI2Y2FBgC/Ilzbo5zblFP7ICIiHRNha+J2cGjHI0VsCEyCuhby0l01R133MGiRYvYv38/b731FuC8jiQDTDrlazqwGpgErANSnYV5ecxsYAMwDngV2Jb88yagnUT5+ksz23T5aJiIiGTe8ePHmR/6lAY3iLXhMbgBUrwumzlzJvfeey9bt27lroB+AlIyK51rOxYAR4ESYE3ydqoxZ4ERwF6gwDm3IvmR5d8Arzjn2s1sL/Cuc+4DM1ttZpXOuboe2RMREUnLhQsXePXVV7nkgqxsv40ofq8jeWLBggU0NzfjPt5NY3wQJ+Kp/nkT6XnplK9mIAf4hMQRriMdjCkCdgJLkrcB/g543zn3KwDn3CtXPWYzMBr4XPkys6XAUoDy8nJqa2vT3JWuaWho6NXtS9doXrKP5qRz5r+wJ+X9a749Ma2xqcalku68XH4vjcfjvPPOO7S1tbEqPJ42ginHeaHquZrr7kv3+5BKOt/XyZMn8/6uw8wLfcpv2idw0eV0mCWVjsd9/rm7sx/SM7LpPSyd8rUVWOKce83M5gI/62DMBBJHxiYBh8zsK8Cxy8ULwMx+DPyExN/KauBH127IObcMWAZQXV3tqqqqOrM/XZKJ55DO07xkH81JZ6QuX6m/h9eP7cz3+vqxHW/v7bff5rPPPuPJJ5/kxZePdut5r5d6n7ujp/Ok2t574bE8lrOXB0KHqGn/ArFeOBKo1052yJZ5SOecrw1AoZltBDY7506mGPMGsCA59iXnXDvwHeDLyZ9sXGdmd5E4Evb3JMraO865Ez2yFyIiclO7d+9m48aNzJw5k7vuusvrOFnjgsvl/fBohvtauSd4FJ2AL73tpke+nHMOeObybTObB6wCxl4+X8s51wYsvuZxczvY5JwupxURkS45c+YMv/nNb6ioqGDhwoVex8k6x+OFbI+UMjVYz5n4YPbFSryOJP1YVxZZ3QJMAVIdARMRkSwTJMarr75KKBTiS1/6En7/wDzB/mZ2RMs4FitgZvAYxb6LXseRfqzT5cs51+Kc2+2c09LAIiJZz3FvsJbGxkZ+//d/n/z8fK8DZTHj/fBoWlyI+4OHCRH1OpD0U7q8kIhIP3abv4HRgXM88MADWXOycTYLE2BNeAx5FuWeYB06/0t6g8qXiEg/NdTamBU8yslYPvfee6/XcfqMBjeYbdEyRgfOMc6fPcsTSP+h8iUi0g8ZceaGjhDH+CA8GrOBtYJ9d+2OjqA+NoTZwaPkW5vXcaSfUfkSEemHpgTqKfa18GG4kkuEvI7T5ziMDyJjiGPMDR3B0AW4peeofImI9DMlvgtMCtRzMDqc2vgwr+P0WS0uxIfhSkp8LUwJ1HsdR/oRlS8RkX6kra2NucEjXHQ5bIiM8jpOn1cbH8bB6HAmBeq51XfB6zjST6h8iYj0I8uXL2ewhXk/MnrAXjC7p22IjOKiy+G+4BECxLyOI/2AypeISD+xf/9+du3axa5oKWfiQ7yO029E8fNBpIohFqY6eNzrONIPqHyJiPQDbW1t/O53v6OkpISd0VKv4/Q7n8Xz+SRWwh2BM/r4UbpN5UtEpB945513aGlp4YknniCut/ZesS1Szvl4DnOCtfj18aN0g16hIiJ93OHDh9mxYwf33HMPZWVlXsfpt6L4WR+pZKivnenBE17HkT5M5UtEpA+LRCL89re/paioiPnz53sdp987FR/K3mgxE/yfUaKLb0sXqXyJiPRhW7dupbm5mccff5xAIOB1nAFhS2QkF12Ie4O1+LX4qnSBXqkiIh2oeq7muvtqf/CIB0lSq62tZf/+/cyePZuKiooe2Wa27zN4nzGKnw8jlTycc5CpgZNsiY7ssW2n2jdIf/+8/t5IenTkS0SkD4pEIvzmN78hPz+fBx54wOs4A87JeAH7o0VMDJximF3yOo70MSpfIiJ90AcffMC5c+e4++67CQaDXscZkLZERtJOgHtCdRjO6zjSh6h8iYj0MWfOnGH9+vVMmjSJ0lKt6eWVMAE2RSoo9rUw3n/G6zjSh6h8iYj0Ic45ampqCIVCLFy40Os4A96nsWGciOUzPXiCPMJex5E+QuVLRKQP2bVrF3V1dTz44IMMHjzY6ziCsSFSiY84M4PHvA4jfYTKl4hIH3Hp0iXeffddRo4cybRp07yOI0nnXS67oqWMCZyjzNfsdRzpA1S+RET6iJUrV9La2sqjjz6KmXkdR67ycXQETfFc7g4e1dpfclMqXyIifcDRo0fZvn07d999N7feeqvXceQacXx8FBnFUF87kwMnvY4jWU7lS0Qky8ViMWpqaigoKGDevHlex5EOnIoP5WB0OHcGTlNgrV7HkSym8iUikuU2b97MZ599xhe/+EVCoZDXceQGtkRGEsXHrOAx0Npf0gGVLxGRLHbx4kXWrFnDuHHjGD9+vNdx5CbaCLI9Uka5/zyjfE1ex5EspfIlIpLFVq5cSSQS4Ytf/KJOsu8j9sVKaIznMTN4DD8xr+NIFlL5EhHJUseOHWPnzp3cfffdDB8+3Os4kiaHsSEyinxfmLsCp7yOI1lI5UtEJAvF43GWL19Ofn4+c+fO9TqOdNLpeD6Ho8O4K3CKIdbudRzJMipfIiJZaNu2bdTX17Nw4UKdZN9HbYmMxGFa+V6uc9PyZQkvmtk2M/teB2NyzOwNM9tuZs8m7xttZqvMbIuZ/VXyvmHJ+3aY2WM9uysiIv1DDlHee+89KisrmThxotdxpIsuEWJHtJRKfxMHDx70Oo5kkXSOfM0GcoFqYKmZlaUY8ySwG5gPPG9mOcB/AP69c64aeMTMhgLPAq8AXwd+2P34IiL9z9TgCdra2li0aJFOsu/jPoneSnM8h7fffptoNOp1HMkS6ZSv6cBqYBKwDph6gzGzgQ3AOOBVYFvyz5uA9qvGTQAazUxXhRURucowu8R4/xlmzJihlez7gTg+NkRG0djYyKZNm7yOI1kikMaYAuAoUAKsSd5ONeYsMALYCxQ451YkP7L8G+AV51y7mV0edw44mXxcy9UbMrOlwFKA8vJyamtru7Bb6WtoaOjV7UvXaF6yj+YkIdV70vwX9qT9+Krnam7wp46ZwaOECfB/rA3z79Z+fuyab1//EWS689Kd99JUmVNl6Q3p5r7x97Xz2+tJJ+MFHI0VUPPuKv7oN5/RRhDo3PewO7l7ep9T/X3P1N+H7sim97B0ylczkAN8QuLI1ZEOxhQBO4ElydsAfwe875z71VXjxpI4gvbdq8Zd4ZxbBiwDqK6udlVVVWnuStdl4jmk8zQv2Wfgzcn1/8ik/h6kX75upNLXRKn/Ih+GRxFO8fbc0ff/+vt7P3dvfh9683kylftamyMV/F7OHqYFT/BhpKrTWdJ/7aU7992RiefoHdmSM52PHbcCs5xzh4C5wPYOxkwAzpD4ePKQmX0FOHZV8bo8rjz5vPnOuZbrtiQiMgD5iTMjeIxz8TwOxIq9jiM97LzLZW+0hNv9Zxlml7yOIx5Lp3xtAArNbCOw2TmX6nLtbwALkmNfcs61A98Bvmxm65K/7iJxROtPSXx8+f2e2AERkf5gQuA0+b4wGyMVOHSSfX+0M1pKOwFmBY+i6z4ObDf92NE554BnLt82s3nAKmCsc64uOaYNWHzN4zpaFXBBl9OKiPRDeYSZHKinLlZIfXyo13Gkl4QJsC1Szj2hOqp857yOIx7qyiKrW4ApJE6YFxGRbpoWPIEPx+bISK+jSC87ECuiMZ5HdfA4kUjE6zjikU6XL+dci3Nut3NOf2tERLppuLVwm7+BT6IlXHC5XseRXuYwNkYqyPeF+eijj7yOIx7R5YVERDzjmBU8RhsBdkZTrV8t/dGp+FDqYoWsW7eOCxcueB1HPKDyJSLikSrfOW71X2RrpJwIfq/jSAZtilQQj8d57733vI4iHlD5EhHxgJ841cHjNMbzOBQr8jqOZNhFl8PMmTPZsWMH9fX1XseRDFP5EhHxwB3JpSU2aWmJAWvu3LkMGjSId999l8TCAjJQqHyJiGRYLhEmB05xNFagpSUGsNzcXObPn09tbS379+/3Oo5kkMqXiEiGTQ2eJECczZEKr6OIx6ZPn05RURErVqwgFot5HUcyROVLRCSDCq2V2/1n2Bcr5ryWlhjwfD4fCxcupLGxkc2bN3sdRzJE5UtEJINmBI8Rwc+OiJaWkIRx48YxduxY1q5dSw5Rr+NIBqh8iYhkSLmvmZH+8+yIJK7xJwJgZixcuJD29namBHXxmIFA5UtEJAMMx4zgMc7Hc9gXK/E6jmSZkpISpk2bxhf8ZxhqbV7HkV6m8iUikgG3+89wi6+NLZGRxPXWKyncf//9RDFmBI97HUV6md4BRER6WZAYU4MnORUbQl280Os4kqUGDx7Mrmgpo/xNjPCd9zqO9CKVLxGRXnZXoJ48i7IpUgFaUFVu4JPorVyMh5JHv7Twan+lMz5FRHrRYAszMXCaw9FhNLjBXsfpc6qeq8nq7XVH6iw+tkbLmRc6wlh/A4d16al+SUe+RER60bRA4vydrdFyj5NIX/FpbBhn4oOYHjyBHy282h+pfImI9JLh1sK4QCN7orfS4nK8jiN9hrE5UsFgi3Bn4LTXYaQXqHyJiPQKx8zgMVpdgI+jpV6HkT7mdDyfulghdwVOkUfE6zjSw1S+RER6wShfEyP8F9keKSOC3+s40gdtiYzEj2Nq8ITXUaSHqXyJiPQwH3Gqg8dpiudyIFbsdRzpo867XPbGirnNf5bTp/XxY3+i8iUi0sPG+89Q4Gtnc2QkTktLSDfsSB45XbFihddRpAepfImI9KAQUaYET3Iyls/xeIHXcaSPCxNgZ7SUw4cPc+jQIa/jSA9R+RIR6UGTA/XkENOCqtJj9kZLuOWWW1ixYgXxeNzrONIDVL5ERHrIEGvnjsBnHIwVcc4N8jqO9BNxfDz44IN89tlnbN++3es40gNUvkREekh14DhxjG2RMq+jSD9zxx13UFFRwerVq2lvb/c6jnSTypeISA8o8V1gdOAcu6MjaCXkdRzpZ8yMhQsX0tLSwvr1672OI92k8iUi0m2OGcHjtLggu6O3eh1G+qmRI0dy55138tFHH3H+/Hmv40g3qHyJiHTTaP85SnwtbIuUE9WCqtKLFixYgHOO9957z+so0g0qXyIi3RCNRpkeOE5DPI/DseFex5F+rrCwkFmzZrFz505OnjzpdRzpIpUvEZFu2LBhA/m+MJsjFVpQVTLivvvuY9CgQbz77rs457yOI11w0/JlCS+a2TYz+14HY3LM7A0z225mzybv85vZ82a27qpx3zGzj81snZkt77ndEBHJvIsXL/LBBx9wNFZAfXyo13FkgMjNzWX+/PnU1dWxb98+r+NIF6Rz5Gs2kAtUA0vNLNXPUD8J7AbmA8+bWQ7gB9YDgavGDQH+xDk3xzm3qDvBRUS8tmbNGqLRKJsjFV5HkQFm+vTpFBcXs2LFCqLRqNdxpJPSKV/TgdXAJGAdMPUGY2YDG4Bxzrmwc67mmnFDgL80s01m9lddjy0i4q3Tp0+zbds2qqurOe9yvY4jA4zP52PhwoWcO3eOzZs3ex1HOilw8yEUAEeBEmBN8naqMWeBEcDeDsaQ/LN3nXMfmNlqM6t0ztVdPcDMlgJLAcrLy6mtrU0jYtc1NDT06valazQv2aevzsn8F/Zcd9+ab0/s8vZqa2txzrFixQqCwSBjxoyBtZm55l6q98PE/l2/j9eqeu7a/wt3T09vT24s1dw/+A/7eSg0lN++s4pvvHmadoJpP7anpfr70J3XWW/IpvewdMpXM5ADfELiCNeRDsYUATuBJcnb13HOvXLVzc3AaKDumjHLgGUA1dXVrqqqKo2I3ZOJ55DO07xkn745J9cXk/T3I/VjDxw4QH19PQ8//DDjx48HMlO+Uue+efGSvq+jud8UqeD3cvYwJVjPxsioTjy2O9L7O5eN7xfZkimdjx23ArOcc4eAuUCqC0ttBSYAZ0h8PJnyncjMfmxmd5mZj8Q5ZAe7lFpExCOxWIwVK1YwbNgwZsyY4XUcGeCaXR77Y8V8wf8ZBdYkk3OdAAAcPElEQVTqdRxJUzrlawNQaGYbgc3OuVQLi7wBLEiOfck519GFp/4O+HsSZe0d59yJLmQWEfHM1q1bOXv2LAsXLsTv14Kq4r3tkTIi+JkRPO51FEnTTT92dIlFRJ65fNvM5gGrgLGXz9dyzrUBizt4/Oyrvt4PzOlmZhERT4SIsmbNGkaPHs3tt9/udRwRANoJsitayozgccp8zZyMd3TatWSLriyyugWYAmhpXREZUCYH62ltbWXhwoWYaUFVyR6fREs4H89hZvAYhhZezXadLl/OuRbn3G7nXKQ3AomIZKMCa2WC/zOmTZvGiBEjvI4j8jlxfGyOjOQWXxtf8J/xOo7chC4vJCKShhnB40Tx8cADD3gdRSSlo/FCTsbymRI8QQgtvJrNVL5ERG6i3NdMhb+ZHdFSBg8e7HUckQ4YmyIVhIgxNagzg7KZypeIyA0YcWYGj9Ecz2FvtMTrOCI3dM4NurL0RKGWnshaKl8iIjdwh/8Mhb42NkUqiOstU/qAy0tPzAweA518n5X0TiIi0oEcIkwJnuR4bCjH9eP70ke0E2RHpIxy/3lG+lJecEY8pvIlItKBacGTBImxKVIBaGkJ6Tv2xoppiucyM3iMWCzmdRy5hsqXiEgKp0+f5nb/GfbFSmh2eV7HEekUh49NkQoKfO1s3LjR6zhyDZUvEZFrOOdYvnw5Yfxsj5R5HUekS07ECzgWK2Dt2rVcuHDB6zhyFZUvEZFr7Nmzh7q6OrZFRhK++VXYRLLWxkgFsViMlStXeh1FrqLyJSJylXA4zLvvvktpaSkHYkVexxHplgsul7vvvptdu3Zx9OhRr+NIksqXiMhV3n//fS5cuMCiRYtwOsle+oH77ruPoUOHsnz5cuLxuNdxBJUvEZErGhoa+Oijj5g8eTIVFRVexxHpEaFQiIULF3Lq1Cm2bdvmdRxB5UtEBPjXk+yDwSAPPvig13FEetSECROoqqrivffe49KlS17HGfB0JqmIZJWq52rSGlf7g0d69Hn379/P4cOHefjhhxkyZEiPbrunpPu9Ebna5b83hTaYJ3Ja+db3X+IXf/2/d+qxmZbqeTvzmk/1+DXfntitTD1JR75EZMCLRCK88847FBcXM2PGDK/jiPSKJpfH3lgJ4/1nqK+v9zrOgKbyJSID3vr162lqamLRokX4/X6v44j0mh2RMtoIUFNTg3O67qNXVL5EZEBrbGxk3bp1TJw4kdGjR3sdR6RXhQmwOTKSEydO6OR7D6l8iciA5Zzjrbfewu/38/DDD3sdRyQjDseGU1lZycqVK2lpafE6zoCk8iUiA9Ynn3zC4cOHeeCBB8jPz/c6jkiGGI888gjhcFgr33tE5UtEBqT29nbefvttRowYoZPsZcApLi7m7rvvZseOHdTV1XkdZ8BR+RKRAWn16tVcvHiRRx99FJ9Pb4Uy8MybN4+CggJqamqIxWJexxlQ9I4jIgNOfX09mzZtorq6mvLycq/jiHgiGAyyaNEizpw5w4YNG7yOM6CofInIAOOoqalh0KBBLFiwwOswIp4aP34848ePZ+3atTQ1NXkdZ8BQ+RKRAWW8/ywnTpxg4cKF5Obmeh1HxHOLFi0CYPny5Vr7K0NUvkRkwMgjzPTgcaqqqrjrrru8jiOSFQoKCpg/fz4HDhxg7969XscZEFS+RGTAuDt0FD9xHnvsMczM6zgiWWP27NmUlpby1ltv0dra6nWcfk/lS0QGhErfOSr9TWyPljFs2DCv44hkFZ/Px+OPP86lS5dYsWKF13H6PZUvEen3QkSZHTpKQ3wQe6IjvI4jkpVGjBjBPffcw/bt2zly5IjXcfq1m5YvS3jRzLaZ2fc6GJNjZm+Y2XYzezZ5n9/MnjezdVeNG2Zmq8xsh5k91nO7ISLSsergcXKJsD5ciUMfN4p0ZN68edxyyy387ne/IxKJeB2n30rnyNdsIBeoBpaaWVmKMU8Cu4H5wPNmlgP4gfVA4KpxzwKvAF8Hftj12CIi6RnhO8/4wFn2RG+lwQ32Oo5IVgsGgzz22GM0Njaydu1ar+P0W+mUr+nAamASsA6YeoMxs4ENwDjnXNg5V9PBuAlAo5npnVBEek0kEuHeYB3n4zlsj6b6f6OIXGv06NFMmTKFDz/8kFvsktdx+qV0ylcBcBYoAdYkb3c0xgF7Oxhz9bhzwMkbjBMR6ba1a9cy1NfOh5FKYvi9jiPSZyxcuJBBgwYxJ1SLEfc6Tr8TuPkQmoEc4BMSR65SnYXXDBQBO4ElydsdbWssiSNo3001zsyWAksBysvLqa2tTSNi1zU0NPTq9qVrNC/ZJ9vm5GbvDWfOnOHDDz/kQLSI+vjQz/1Z1XPXHpSHNd+e2CPPK9IbevrvXTrbmzFjBi1r1jApcIqdXThy3J3XWSrd/R5k03tYOuVrK7DEOfeamc0FftbBmAkkjoxNAg7dYFvlwAEg3znXcu0A59wyYBlAdXW1q6qqSiNi92TiOaTzNC/ZJzNzsietUTfKEolEqKmpIT8/n02nR3Zje9dnSXecSE/q6b936byWq6qq+MeVu5gSqOdYrJBGN6jLz9eZ501I97WX/uOHDx+eNf+upPOx4wag0Mw2ApudcydTjHkDWJAc+5Jzrr2DbS0D/pRESft+5+OKiNzc6tWrOXv2LI8//jiRtP6PKSKpbIyMoo0Ac0JH8Onjxx5z03cll7jQ0zOXb5vZPGAVMNY5V5cc0wYs7uDxs6/6uoFESRMR6RVHjx7lo48+Yvr06YwdOxbY53UkkT6rnQAfhit5MOcQkwL17IiWex2pX+jKIqtbgCkkTpgXEcka4XCYN998k8LCQh566CGv44j0C8fihRyKDmdyoJ7hdt3ZQtIFnS5fzrkW59xu55xWXxORrLJq1SoaGxt54oknyMnJ8TqOSL+xMVJBK0HuC9Xq48ceoMsLiUi/cOTIETZt2sTMmTOz5qRakf4iTID14Spu8bUyJaAPvrpL5UtE+ry2tjbefPNNhg0bxoIFOq1UpDeciBdwIFrEXYFT3Oq74HWcPk3lS0T6NOccNTU1nD9/nsWLFxMKhbyOJNJvbYxUcNHlMDd4hBBRr+P0WSpfItKn7dq1i927dzN//nxGjkxvTS8R6ZooftaGRzPIwtwdPEriwjbSWSpfItJnnTt3jrfeeotRo0YxZ84cr+OIDAhn3RC2R8sYE2hkrL/R6zh9ksqXiPRJ8XicX//615gZixcvxufT25lIpnwcLeVUbAizg3UMsY7WVZeO6N1KRPqktWvXcvz4cR599FEKCwu9jiMyoDiM9yOjcRjzgp9i+vixU1S+RKTPKfFd4IMPPmDy5MnceeedXscRGZBaXA4fRSop8bcwWctPdIrKl4j0KSGizAseobCwkEWLFnkdR2RAOxIbdmX1+1Lfea/j9BkqXyLShzjmho6QZxGWLFmiVexFssBHkVGcd7nMDX1KHmGv4/QJKl8i0mfcFThFhb+ZTZEKyst1gV+RbBDFz+rwWELEmRfS+V/pUPkSkT5hhO880wIn+DQ6jH2xYq/jiMhVmlweH0ZGUeq/yNTACa/jZD2VLxHJenlEmBf6lPMul/WRSsC8jiQi1zgcK2J/tIjJwVOM9DV5HSerqXyJSFYzHPNCnxIizurwWKL4vY4kIh3YGBlFQzyP+0JHGKz1vzoU8DqAiPQvVc/VpD229geP3HTM1MAJSv0XeD9cRZPLu+HzpLM9kb6oM68rL58jho814bE8lvMJ80Ofsrx9PPGrjvPodZugI18ikrUqfeeYHDzF/mgRh2NFXscRkTScd7msC4+mxNei6z92QOVLRLLSLXaJ+0JH+Cw+mA2RUV7HEZFOqIvfwo5IKbcHznKH/zOv42QdlS8RyTo5RFgQOkQYP++1j/3cxxYi0jdsj5ZRFytkZvCYFmC9ht7RRCSr+IjzQOgweRZhVfs4Wgl5HUlEusR4PzyaZpfL/NBhXYD7KipfIpJVZgWPMcJ/kfWRKhrcYK/jiEg3RPGzKjwOAx4MHSRAzOtIWUHlS0SyxubNm/lC4Ay7IiP4NDbc6zgi0gMuuFxWh8dSYG3MDR1BJ+CrfIlIljh06BBvv/02x2IFbIvq0kEi/Ul9fCibIhVU+puYETzudRzPqXyJiOfq6+t57bXXKC4uZm14DE4r2Iv0O3tjJeyJlnBn4DQT/Ke9juMplS8R8VRTUxOvvPIKeXl5fO1rXyOiFexF+iljc6SC2uRPQFb6znkdyDMqXyLimdbWVn7+858TjUb52te+Rn5+vteRRKQXOYz3w2P4LD6YuaFPKfFd8DqSJ1S+RMQTfuL84he/oKmpia9+9asUFxd7HUlEMiCGj1XhcbS4EAtChzh79qzXkTJO5UtEMs5wzA0d4dixYyxevJjKykqvI4lIBrUT5N3w7TiMl19+mQsXBtYRMJUvEckwx73BWqr853j44YeZOHGi14FExAMXXQ4r2m/j0qVL/NM//RMtLS1eR8oYlS8RySDH7OBRbgs0sD1SyuzZs70OJCIeanCD+YM/+AOampr453/+Z1pbW72OlBE3LV+W8KKZbTOz73UwJsfM3jCz7Wb2bPK+YWa2ysx2mNljyft+aGabzWydmf20Z3dFRLKbY0bgOHcEzvBxZAQ7omVeBxKRLFBZWclXv/pVzp49y8svv0x7e/+/DFE6R75mA7lANbDUzFK9Yz4J7AbmA8+bWQ7wLPAK8HXgh8lxQ4AvOefmOOe+2c3sItKHTA2c5M7gaT6JlrAlWg5ay0tEksaOHcuXvvQl6uvreeWVV/D388sQpVO+pgOrgUnAOmDqDcbMBjYA4666bwLQaGaDSZSvvzWzrWa2tPvxRaQvmBSoZ0qwnv3RIjZGKlDxEpFrjR8/nsWLF3Ps2DEWhA7jJ+51pF6TTvkqAM4CJcCa5O2Oxjhgb/L25fvOASeTt7cCf06ipH3HzHK7F19EsptjcuAk04MnOBQdxkeRSlS8RKQjd955J48//jhlvvP9+kLcgTTGNAM5wCckjmYd6WBMEbATWJK83QyMJXG07LtAs3Pu/738ADPbD4wAaq/eUPKI2FKA8vJyams/98c9rqGhoVe3L12jeck+qeZk/gt7bvAIR3XgBHcFT3EoOpx1karrLhvU06/vqudquvzYdLP09nuSSCZ057XS0659TRUWFvJBZDRzgkdYGDrAyvBthAl0mHnNt9P7iels+nclnfK1FVjinHvNzOYCP+tgzAQSR8YmAYeS95UDB4B851yLmf0LiXPBLgJjgPprN+ScWwYsA6iurnZVVVWd3KXOy8RzSOdpXrLP9XPSUflK/FTjHYEz7I0WsyEyilRHvFLP8Y0KXe9JN0s2ZRbpD1K9pg7H9hB1PuaFPuXhnAO8234b7QTTfnyq1+Tw4cOz5t+VdD523AAUmtlGYLNz7mSKMW8AC5JjX3LOtZMoUH9KopB9Pznu/wFqkuN+lBwnIv2I4bgvWHvlpxo7Kl4iIjdSF7+F98JjKbRWFuXsJ4+w15F6zE2PfDnnHPDM5dtmNg9YBYx1ztUlx7QBi695XAOJQnb1feuBWd2PLSLZyEeceaFPqfI3sTVSxq5oKSpeItJVx+OFrAjfxoOhQ/ybnP28E76diy7H61jd1pVFVrcAU0icRC8iAkAOUR4OHaDK38TGcAW7omWoeIlId52KD+Wd9tvJsSiP5uylyC56HanbOl2+nHMtzrndzrlIbwQSkb4n39p4JGcvRb4W1oTH8EnsVq8jiUg/csYN4XftdxBxfhbl7KfSd87rSN2iywuJSLeU+C7waM4+cizGO+HbORIb5nUkEemHzrtcftf+BRrdIO4PHWaC/xSJFa76nnR+2lFEJKXR/kbuCx7hoguxInwbF5yW7hOR3tNOkLfbxzM3dIRZoeMMjbYTj8fx+frWsaS+lVZEskI8Hue9995jfuhTzsQHU9N+h4qXiGREDB+rw2P4OHIrdwTO8POf/5yWlhavY3WKypeIdMqlS5d4+eWX+eCDDzgQLeKd8O206yC6iGSUsSVawQfhKo4ePcqyZcs4fvy416HSpvIlImk7ceIEP/nJT6irq+Oxxx5jfaSKuN5GRMQjh2JFfOtb38Ln8/HSSy+xZcsWEitkZTf9d1VEbso5x/79+9m8eTP5+fn80R/9EWVlZfCr7LlEiYgMTKWlpSxdupRf//rX1NTUcPz4cfz4iOH3OlqHVL5E5IZaW1upqalhz549jBs3jieffJK8vDyvY4mIXJGXl8dTTz3F2rVref/993k8J5f3w6NpcIO9jpaSypeIdOjw4cO8+eabtLS0MG3aNB599FHMtHCqiGQfn8/H/fffT2VlJS/80y95NGcf26Jl7I6OwGXZgs8qXyJynUgkwsqVK9m0aRNFRUU89dRTtLe3q3iJSNYbM2YM/6ttIvcE66gOnmCkr5kPIqO9jvU5Kl8i8jnHjx/nzTff5OzZs8yaNYsFCxYQDAapra31OpqISFrCBFgTGcOxeCOzg3U8kbOH2tphVFVVeR0NUPkSkaS2tjZWrVrFli1byM/P5+tf/zpjx471OpaISBcZh2PDOR0fwn3BIwSDQa8DXaHyJTLAOefYs2cP77zzDi0tLcyaNYv777+fnJwcr6OJiHTbRZfD8vB4/qK83OsoV6h8iQxgDQ0NLF++nMOHD1NWVsZTTz2VWEJCRKRfya7zVVW+pN+peu76tadqf/BIrz9HR1I9d3cy9sT+tbS0sHbtWrZu3UpbDLZFKth3uIS/+dvtwPaU25v/wh5gT6eepz/qzNyLSO/ri69JlS+RASQSibBhwwbWr19POBxm2rRpfHddlDay51wIEZH+TuVLZACIxWLs3LmTtWvXcv78ecaPH8+CBQsoLi7mj9f1vf81ioj0ZSpfIv1YNBpl27ZtrF+/nvPnz1NeXs7ixYuz5setRUQGIpUvkX4oHA6zdetWPvzwQy5evEhFRQWPPfYYY8eO1UKpIiIeU/kS6UcaGxvZvHkz27dvp729ndGjR7NkyRIqKytVukREsoTKl0if5zh06BCbNm3i4MGD+Hw+JkyYwKxZsxg5cqTX4URE5BoqXyJ91BBrZ5y/gbH+Bl5+eSuDBw9m7ty5VFdXk5+f73U8ERHpgMqXSB/S3t7Obf4zjPM3MMJ/EeegPp7P/7b43zBhwgQCAb2kRUSynd6pRbLcpUuX2L9/P3v37uXTTz9lTihGczyHrZFyDseG0eJy+MmkSV7HFBGRNKl8iWQZ5xyNjY0cOnSI/fv3U1tbi3OOgoICZsyYwf+9pokzbjDZdrkMERFJj8qXSBYIEWXv3r0cOnSIw4cP09zcDEBRURFz5szhjjvuYMSIEZgZz67WoqgiIn2ZypeIB/IIM8J/kRLfRW71XWCYtfKrX+0gFAoxZswY7r33XsaNG8ctt9zidVQREelhKl8ivay1tZX6+npOnjx55fev5jUBEHE+zsQHsyNWyo/+7cOMHDkSv9/vcWIREelNKl8iPcRHnAJro9DXyi1Xfm/lv/7XLVfGFBYWUlZWxjufDeF0PJ9Gl4fDB0BlZaVX0UVEJINUvkQ6IRwO09zczEhfE/nWzlBfO0OtjaHWzhBrx5c8Bz7u4LzLpSE+iMUPzaG0tJSysjLy8vIA+POtOm9LRGSgumn5ssQ1SZYB04HfOef+U4oxOcAvgSrgBefcT8xsGPAaMBz4j86535rZGODnQB6w1Dm3ucf2RKQb4vE4ly5doqWlhYsXL3LhwoXP/WpubqapqYnW1lYAHspJPC7sfFdK1hE3jHMul6Z4Hs0ul3jyiNZLc+Z4tVsiIpKF0jnyNRvIBaqBk2b2gnPu5DVjngR2A98APjWzl4BngVeAjcDrwG+BvwD+I4ny9dfAF7u/CyIJsViMcDjMYGsnRIygxa78vmXLFlpbW6/8amtro7W1lZaWFlpaWq6Uqmvl5uaSn59PQUEBZWVlFBQUUFhYyLd+sYfz8VzaCKAlH0REpDPSKV/TgdXAJGAdMBW4tnxNB94iUdQ2AOOS9/1fJEpbo5kNBiYDfwJ8BRjWA/mlA865675O5/fLv669nepXPB6/7uurf+/oVywWu+73a381Njayd+9eotEo0WiUWCxGJBIhEokQjUavfB2JRAiHw4TDYeLxOABfzr3++1FTcwSAQCBAXl7elV/FxcVUVVUxaNAgBg8ezODBg8nPzyc/P58hQ4YQDAZTfn8/e/lod6dIREQGqHTKVwFwFCgB1iRvpxpzFhgB7E3evnzfORJlrYDEIQIfUAvEupW8B+zfv5/XX3+dxCerPePq0tNTj011f6py1Z/4/X6CwSCBQIBAIHDl9uVfeXl5V74OhUKEQqErX//H3+4j7PxE8F/5fd1fPkxeXp4uvyMiIp6zm/3DbWZ/DJwGtgNfAnY552quGfPfSBz52gN8D/hb4K+AvwH2Af+LxEeTr5P46PFj4APn3HUnw5jZUmBp8uZ4YH8X9y1dRSRKomQXzUv20ZxkJ81L9tGcZKdMzEulc674ZoPSOQywFVjinHvNzOYCP+tgzAQSR8YmAYeS95UDB4B851yLmR1KPmcZcD7VkznnlpE4wT8jzGyLc646U88n6dG8ZB/NSXbSvGQfzUl2yqZ58aUxZgNQaGYbgc0pTrYHeANYkBz7knOunUSB+lMShez7yXH/A/ifwK+B/9yt5CIiIiJ90E2PfLnE55LPXL5tZvOAVcBY51xdckwbsPiaxzWQKGRX3/cpcHf3Y4uIiIj0Tekc+brWFmAK1//EY1+VsY84pVM0L9lHc5KdNC/ZR3OSnbJmXm56wr2IiIiI9JyuHPkSERERkS4a0OXLEl40s21m9j2v8wxUZpZjZq+Y2Toz+03ytuYlS5jZ75tZvV4v2cHMnjOzzWb2hpkFNSfeS75nvZx8D1tuZrdqXrxjZk+Z2Vkzy031vpWcrzfMbLuZPetFxgFdvvj8pZOWmlmZx3kGqi8BO5Lrvh0EnkPzkhXMbATwFHAMvV48Z2YjgTucczOAbcB30ZxkgweBY8n3sI0k1qrUvHjnAv+6Rmiq963Ll0ScDzyfvD51Rg308pXq0kmSeTtJXAcUIAo4NC/Z4r8DfwbE0eslGzwINJvZuyTWUWxBc5INaoH/08z2k7hm8Tk0L55xzv0OiCRvpnrfunzf1ZdEzKiBXr4uXwLpRpdOkl7mnPvYOXfczH6PxBzE0Lx4zsz+LfDu5SVl0OslG4wAKpxzC0ksVH0LmpNsUA8sds6NB94EhqN5yRap3rcu3+f410siZtRAL1/NQA7wKZCfvC0eMLOvA3OAf4fmJVs8AXzTzNaQuILFn6F58dpFYG3y68u/a0689+ckChgkLqf3n9G8ZItU/540k7jU0E48mp+BXr62ArOcc4eAuSSuXykZZmaVJC5h9d3kor6alyzgnHvUOTffOTcf+AR4FM2L1z4CZiW/ngaE0Zxki3uSv98N/Cc0L9ki1b8nly+JeIZ/vSRiRg308pXOpZOk930LuCv5k0LrgNvRvGQjvV485pzbCpw0s49InKfy39CcZIP/ATyZnJfHgL9F85ItUr1vpbokYkZpkVURERGRDBroR75EREREMkrlS0RERCSDVL5EREREMkjlS0RERCSDVL5EREREMkjlS0RERCSDVL5EREREMkjlS0RERCSD/n+ul96usEXP6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "xs = np.arange(101)\n",
    "rv = stats.norm(p_mean, np.sqrt(p_var))\n",
    "ax.plot(xs, rv.pdf(xs), color='gray')\n",
    "ax.hist(scores, bins=100, range=(0, 100), density=True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.539373Z",
     "start_time": "2018-08-14T18:04:47.535604Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([49, 60, 72, 86, 52, 61, 77, 91, 80, 56, 69, 67, 90, 56, 75, 79, 60,\n",
       "       79, 68, 81])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "n = 20\n",
    "sample = np.random.choice(scores, n)\n",
    "sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.551920Z",
     "start_time": "2018-08-14T18:04:47.540558Z"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(1111)\n",
    "n_samples = 10000\n",
    "samples = np.random.choice(scores, (n_samples, n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 点推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 母平均の点推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.556646Z",
     "start_time": "2018-08-14T18:04:47.553222Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1回目の標本平均: 67.000\n",
      "2回目の標本平均: 72.850\n",
      "3回目の標本平均: 69.200\n",
      "4回目の標本平均: 64.450\n",
      "5回目の標本平均: 72.650\n"
     ]
    }
   ],
   "source": [
    "for i in range(5):\n",
    "    s_mean = np.mean(samples[i])\n",
    "    print(f'{i+1}回目の標本平均: {s_mean:.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.561429Z",
     "start_time": "2018-08-14T18:04:47.558041Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.538"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_means = np.mean(samples, axis=1)\n",
    "np.mean(sample_means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.577846Z",
     "start_time": "2018-08-14T18:04:47.562500Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.543"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(np.random.choice(scores, int(1e6)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.581915Z",
     "start_time": "2018-08-14T18:04:47.579010Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "70.400"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_mean = np.mean(sample)\n",
    "s_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 母分散の点推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.586398Z",
     "start_time": "2018-08-14T18:04:47.583206Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1回目の標本分散: 116.800\n",
      "2回目の標本分散: 162.928\n",
      "3回目の標本分散: 187.060\n",
      "4回目の標本分散: 149.148\n",
      "5回目の標本分散: 111.528\n"
     ]
    }
   ],
   "source": [
    "for i in range(5):\n",
    "    s_var = np.var(samples[i])\n",
    "    print(f'{i+1}回目の標本分散: {s_var:.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.593666Z",
     "start_time": "2018-08-14T18:04:47.587634Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "196.344"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_vars = np.var(samples, axis=1)\n",
    "np.mean(sample_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.598946Z",
     "start_time": "2018-08-14T18:04:47.594834Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "206.678"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_u_vars = np.var(samples, axis=1, ddof=1)\n",
    "np.mean(sample_u_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.619062Z",
     "start_time": "2018-08-14T18:04:47.600152Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "207.083"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.var(np.random.choice(scores, int(1e6)), ddof=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.624323Z",
     "start_time": "2018-08-14T18:04:47.620579Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "158.253"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_var = np.var(sample, ddof=1)\n",
    "u_var"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 区間推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 正規分布の母平均(分散既知)の区間推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.634102Z",
     "start_time": "2018-08-14T18:04:47.626189Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(64.100, 76.700)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = stats.norm()\n",
    "lcl = s_mean - rv.isf(0.025) * np.sqrt(p_var/n)\n",
    "ucl = s_mean - rv.isf(0.975) * np.sqrt(p_var/n)\n",
    "\n",
    "lcl, ucl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:47.815248Z",
     "start_time": "2018-08-14T18:04:47.637157Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAJBCAYAAACJaIt7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAH8RJREFUeJzt3W+IZHt+1/HPLz2VLJuQSmAT0GnNRQM+iBVvOoWuKLHhRNTIBCFoiKvYErkI/rtCoj4SFRIEH9i0iror2hIWlIhyKSWKaXIfBJ2Ybm0tFWIWM1maoGYjHOMutVvp+vlgZtbp2z3T85vb1af/vF6w3KnfnK7zHaZu77t+51TfUmsNAACv7yuGHgAA4LYRUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANHqw7hN87GMfq2+99da6TwPcMV/60pfylV/5lUOPAdwzR0dHn6u1fsNlx609oN56660cHh6u+zTAHfPkyZN48wVct1LKz7/OcS7hAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAECjB0MPAPBB29vbWSwWefz48dCjAFzIDhQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0enDZAaWUr0ryD5L82iT/O8nvT/I3k3x7kn9ea/2La50QeKX5fJ6Dg4P0fZ/xeJyu6zKZTIYeC+BOe50dqN+f5LjW+tuT/GySv5DkI0mmSd4ppfzqNc4HvMJ8Ps9sNkvf90mSvu8zm80yn88Hngzgbrt0ByrJf0zy/rNf/0qSmuQnknxrkp9M8m1JfmEdw/H69vf3hx6BAZycnOT09PTM2nK5zHvvvZejo6OBpvrw3n777axWK6/rRjs7O0OPAPfGpQFVa50nSSnl9yUZJ/k/ST6X5BvzNKzGH/yaUso7Sd5JkocPH+bJkydXNjAXWywWQ4/AAD4YTy+u3+bXxGq1Sq31Vv8ZhuB7LVyfUmu9/KBS/lCSt5P8YJI/keR/JvkPeXp57z/VWv/Fy752Op3Ww8PDq5kWOGN3d/fLl+9eNB6P8+677w4w0dXY3t7OYrHI48ePhx4FuGdKKUe11ullx116D1Qp5ZuSfE+t9Qfq09o6SvJbaq2fSfIdeRpSwAC6rstoNDqzNhqN0nXdQBMB3A+vcw/U9yeZlFJ+8tnjTyX5ulLKTyX5V7VW9z/BQJ5/2s6n8ACu1+vcA/UXk3zwRxX8w/WMA7SaTCaCCeCa+UGaAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADR6MPQAALzcfD7PwcFB+r7PeDxO13WZTCZDjwX3noACuKHm83lms1mWy2WSpO/7zGazJBFRMLBSa13rCabTaT08PFzrOeCu29/fH3qEa3V8fJzVapWtra2hRxnUyclJTk9Pz61vbGxkc3NzgIkYys7OztAj3BullKNa6/Sy49wDBXBDXRRPr1oHro9LeHAL3Ld3n9vb21ksFtnb2xt6lEHt7u6m7/tz6+Px+N69JuCmsQMFcEN1XZfRaHRmbTQapeu6gSYCnrMDBXBDPb9R3Kfw4OYRUAA32GQyEUxwA7mEBwDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjR4MPQBw/8zn8xwcHKTv+4zH43Rdl8lkMvRYAK9NQAHXaj6fZzabZblcJkn6vs9sNksSEQXcGqXWutYTTKfTenh4uNZzcDvt7+8PPQIDODk5yenp6bn1jY2NbG5uJkmOj4+zWq2ytbV13eNxzXZ2doYeAc4opRzVWqeXHeceKOBaXRRPr1oHuIlcwmMw3nneT7u7u+n7/tz6eDz+8mtie3s7i8Uie3t71zwdwOuxAwVcq67rMhqNzqyNRqN0XTfQRADt7EAB1+r5jeI+hQfcZgIKuHaTyUQwAbeaS3gAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANHow9ABw3ebzeQ4ODtL3fcbjcbquy2QyGXosAG4RAcW9Mp/PM5vNslwukyR932c2myWJiALgtZVa61pPMJ1O6+Hh4VrP8Sb29/eHHoEBnJyc5PT09Nz6xsZGNjc3B5iIixwfH2e1WmVra2ut59nZ2Vnr8wO3TynlqNY6vew490Bxr1wUT69aB4CL3NtLeN553k+7u7vp+/7c+ng89pq4Qba3t7NYLLK3tzf0KAAXsgPFvdJ1XUaj0Zm10WiUrusGmgiA2+je7kBxPz2/Udyn8AD4MAQU985kMhFMAHwoLuEBADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAECjB0MPAMBT8/k8BwcH6fs+4/E4XddlMpkMPRZwAQEFcAPM5/PMZrMsl8skSd/3mc1mSSKi4AYqtda1nmA6ndbDw8O1ngPumv39/aFHGNTx8XFWq1W2traGHuXanJyc5PT09Nz6xsZGNjc3B5iIIe3s7Aw9wr1VSjmqtU4vO849UAA3wEXx9Kp1YFgu4cENdN/ffW5vb2exWGRvb2/oUa7N7u5u+r4/tz4ej+/96wFuIjtQADdA13UZjUZn1kajUbquG2gi4FXsQAHcAM9vFPcpPLgdBBTADTGZTAQT3BIu4QEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQ6MFlB5RSvi/J30iymeTXJ/mxJJ999ts7tdbPrG884LrM5/McHByk7/uMx+N0XZfJZDL0WAA30qUBleSXk/zMs19/TZK/X2v9S2ubCLh28/k8s9ksy+UySdL3fWazWZKIKIALlFrr5QeV8n6S353ktyX54SS/kuQXk3xvrfWLr/ra6XRaDw8PP/ykXIv9/f2hR2AAJycnOT09Pbe+sbGRzc3Na5/n+Pg4q9UqW1tb135ukp2dnaFHgMGUUo5qrdPLjnudHagXfS7Jfq31b5dS/nKSP5DkRy44+TtJ3kmShw8f5smTJ42nYSiLxWLoERjARfH0fH2I18RqtUqt1etxIL5nw+WadqBqrYsX1r4ryW++7HKeHSi4+XZ3d9P3/bn18Xicd99999rn2d7ezmKxyOPHj6/93MD99ro7UE2fwiulfF8p5Y8/e/jx/P97o4BbrOu6jEajM2uj0Shd1w00EcDN1vpjDP5Zku8qpTxO8huS/OjVjwRct8lkkkePHmU8Hid5uvP06NEjN5ADvMRr3QNVa91+4eF3r2cUYEiTyUQwAbwmP0gTAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaPRh6AODNzefzHBwcpO/7jMfjdF2XyWQy9FgAd56AgltqPp9nNptluVwmSfq+z2w2SxIRBbBmpda61hNMp9N6eHi41nPcd/v7+0OPwABOTk5yenp6bn1jYyObm5sDTHR1jo+Ps1qtsrW1NfQot8bOzs7QI8CdUEo5qrVOLzvOPVBwS10UT69aB+DquIR3B3jneT/t7u6m7/tz6+Px+Na/Jra3t7NYLLK3tzf0KAAXsgMFt1TXdRmNRmfWRqNRuq4baCKA+8MOFNxSz28U9yk8gOsnoOAWm0wmgglgAC7hAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANDowdADANxU8/k8BwcH6fs+4/E4XddlMpkMPRZwAwgogAvM5/PMZrMsl8skSd/3mc1mSSKigJRa61pPMJ1O6+Hh4VrPAeu0v78/9Aj3zvHxcVarVba2tgab4eTkJKenp+fWNzY2srm5OcBEDGlnZ2foEbgmpZSjWuv0suPcAwVwgYvi6VXrwP3iEh5cwjvP67e9vZ3FYpG9vb3BZtjd3U3f9+fWx+Ox1wRgBwrgIl3XZTQanVkbjUbpum6giYCbxA4UwAWe3yjuU3jARQQUwEtMJhPBBFzIJTwAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGj0YegCAu2o+n+fg4CB932c8Hqfrukwmk6HHAq6AgAJYg/l8ntlsluVymSTp+z6z2SxJRBTcAaXWutYTTKfTenh4uNZzwE22v78/9Ai3zvHxcVarVba2toYe5Y2dnJzk9PT03PrGxkY2NzcHmIgh7ezsDD0Cr6mUclRrnV52nHugANbgonh61Tpwu7iEB2vmnWe77e3tLBaL7O3tDT3KG9vd3U3f9+fWx+Ox1wTcAXagANag67qMRqMza6PRKF3XDTQRcJXsQAGswfMbxX0KD+4mAQWwJpPJRDDBHeUSHgBAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANHow9ADAzTOfz3NwcJC+7zMej9N1XSaTydBjAdwYAgo4Yz6fZzabZblcJkn6vs9sNksSEQXwTKm1rvUE0+m0Hh4ervUcrMf+/v7QIzCAk5OTnJ6enlvf2NjI5ubmtcxwfHyc1WqVra2tazkfL7ezszP0CHCtSilHtdbpZce5Bwo446J4etU6wH3kEh4v5Z3n/bS7u5u+78+tj8fja3tNbG9vZ7FYZG9v71rOB9DKDhRwRtd1GY1GZ9ZGo1G6rhtoIoCbxw4UcMbzG8V9Cg/g5QQUcM5kMhFMAK/gEh4AQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADR6rYAqpXxfKeVzpZSPlKc+VUr596WUv7LuAQEAbpoHr3ncLyf5mWe//niSjySZJvmFUsrfqbX+wjqGg8vM5/McHByk7/uMx+N0XZfJZDL0WADcca+1A1Vr/edJls8efnuSn0jyrUl+Msm3rWc0eLX5fJ7ZbJa+75Mkfd9nNptlPp8PPBkAd93r7kC9aJzks0m+Mcn7zx4Pbn9/f+gRuGYnJyc5PT09s7ZcLvPee+/l6OhooKm4Cm+//XZWq1XTv9c7Oztrmwfgg94koPokX5Xkv+bpbtTPffCAUso7Sd5JkocPH+bJkycfYsTXs1gs1n4ObpYPxtOL614Pt9tqtUqttenv8Tq+zwA8V2qtr3dgKe8n+d15esnue2qtP1BK+bEk3/+qe6Cm02k9PDy8ilnhjN3d3S9fvnvReDzOu+++O8BEXJXt7e0sFos8fvx46FGAe6aUclRrnV523Jv8GIPHSb6ulPJTSX7aDeQMpeu6jEajM2uj0Shd1w00EQD3xWtfwqu1br/w8I9d/SjQ5vmn7XwKD4Dr9ib3QMGNMZlMBBMA185PIgcAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCoAb7dOfTt56K/mKr3j6z09/euiJIHkw9AAA8DKf/nTyzjvJF77w9PHP//zTx0nyiU8MNxcIKOBW2N4eegKG8Phx8sUvnl37wheS7//+5FOfGmYmhvP++0NP8P+5hAfAjfXBeLpsHa6LHSjgVrhJ7zy5Pm+99fSy3Qd90zd5TTAsO1AA3Fg/9EPJRz96du2jH326DkMSUADcWJ/4RPLJTz7dcSrl6T8/+Uk3kDM8l/AAuNE+8QnBxM1jBwoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARg+GHgBYj/l8noODg/R9n/F4nK7rMplMhh4L4E4QUHAHzefzzGazLJfLJEnf95nNZkkiogCuQKm1rvUE0+m0Hh4ervUcvNz+/v7QIzCAk5OTnJ6enlvf2NjI5ubmABO1OT4+zmq1ytbW1tCj3Ak7OztDjwC3RinlqNY6vew490DBHXRRPL1qHYA2LuHdcd553k+7u7vp+/7c+ng8vhWvie3t7SwWi+zt7Q09CsCF7EDBHdR1XUaj0Zm10WiUrusGmgjgbrEDBXfQ8xvFfQoPYD0EFNxRk8lEMAGsiUt4AACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQ6MHQAwDcZPP5PAcHB+n7PuPxOF3XZTKZDD0WMDABBfAS8/k8s9ksy+UySdL3fWazWZKIKLjnSq11rSeYTqf18PBwreeAddvf3x96hHvl+Pg4q9UqW1tbg85xcnKS09PTc+sbGxvZ3NwcYCKGtLOzM/QIXINSylGtdXrZce6BAniJi+LpVevA/eESHrwG7zyv1/b2dhaLRfb29gadY3d3N33fn1sfj8deE3DP2YECeImu6zIajc6sjUajdF030ETATWEHCuAlnt8o7lN4wAcJKIBXmEwmggk4xyU8AIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABo9GHoAIJnP5zk4OEjf9xmPx+m6LpPJZOixAHgJAQUDm8/nmc1mWS6XSZK+7zObzZJERAHcUKXWutYTTKfTenh4uNZz3BX7+/tDj8AATk5Ocnp6em59Y2Mjm5ubA0w0vOPj46xWq2xtbQ09yrXa2dkZegS490opR7XW6WXHuQcKBnZRPL1qHYDhuYR3g3j3eT/t7u6m7/tz6+Px+N6+Jra3t7NYLLK3tzf0KAAXsgMFA+u6LqPR6MzaaDRK13UDTQTAZexAwcCe3yjuU3gAt4eAghtgMpkIJoBbxCU8AIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBo9OBNvqiU8i1JfizJZ58t7dRaP3NlUwF8wHw+z8HBQfq+z3g8Ttd1mUwmQ48F3FNvFFBJvibJ36+1/qUrnAXgQvP5PLPZLMvlMknS931ms1mSiChgEB8moH5PKeV3JvnFJN9ba/3i1Y0FF9vf3x96BK7B22+/ndVq9eW/75OTk5yenp45Zrlc5r333svR0dEAEzKknZ2doUeANw6ozyXZr7X+7VLKX07yB5L8yPPfLKW8k+SdJHn48GGePHnyYeeEJMlisRh6BK7BarVKrfXLf98fjKfnTk9PvSbuIf+fwk1Qaq0f7glK+a4kv/lll/Om02k9PDz8UOcA7pft7e0sFos8fvw4SbK7u5u+788dNx6P8+677173eMAdVko5qrVOLzvujT6FV0r5vlLKH3/28ONJfuZNngfgdXRdl9FodGZtNBql67qBJgLuuzf9MQb/LMl3lVIeJ/kNSX706kYCOGsymeTRo0cZj8dJnu48PXr0yA3kwGDe6B6oWusiyXdf8SwALzWZTAQTcGP4QZoAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANHow9ADA/Tafz3NwcJC+7zMej9N13dAjAVxKQAGDmc/nmc1mWS6XSZK+7zObzfL5z38+GxsbA08H8HICihthf39/6BEYwMnJSU5PT8+sLZfLTKfTjEYjr4t7ZmdnZ+gR4LW5BwoYzAfj6bla6zVPAtDGDhQ3gnee99Pu7m76vj+3/uM//uP52Mc+lr29vQGmAricHShgMF3XZTQanVkbjUb5+q//+oEmAng9AgoYzGQyyaNHjzIej5Mk4/E4jx49yld/9VcPPBnAq7mEBwxqMplkMpkMPQZAEztQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQ6MHQAwDcRfP5PAcHB+n7PuPxOF3XZTKZDD0WcEUEFMAVm8/nmc1mWS6XSZK+7zObzZJERMEdUWqtaz3BdDqth4eHaz0H3FT7+/tDj3ArHR8fZ7VaZWtra+hR3sjJyUlOT0/PrW9sbGRzc3OAiRjSzs7O0CPQoJRyVGudXnace6AArthF8fSqdeD2cQkP1sg7zzezvb2dxWKRvb29oUd5I7u7u+n7/tz6eDz2moA7wg4UwBXrui6j0ejM2mg0Std1A00EXDU7UABX7PmN4j6FB3eXgAJYg8lkIpjgDnMJDwCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGD4YeAICXm8/nOTg4SN/3GY/H6bouk8lk6LHg3hNQADfUfD7PbDbLcrlMkvR9n9lsliQiCgZWaq1rPcF0Oq2Hh4drPQfcdfv7+0OPcK2Oj4+zWq2ytbU19CiDOjk5yenp6bn1jY2NbG5uDjARQ9nZ2Rl6hHujlHJUa51edpx7oABuqIvi6VXrwPVxCQ9ugfv27nN7ezuLxSJ7e3tDjzKo3d3d9H1/bn08Ht+71wTcNHagAG6orusyGo3OrI1Go3RdN9BEwHN2oABuqOc3ivsUHtw8AgrgBptMJoIJbiCX8AAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoNGDoQcAuCrz+TwHBwfp+z7j8Thd12UymQw9FnAHCSjgTpjP55nNZlkul0mSvu8zm82SREQBV67UWtd6gul0Wg8PD9d6DnjR/v7+0CPwIR0fH2e1WmVra+u1v+bk5CSnp6fn1jc2NrK5uXmV43EL7OzsDD0Ct1Qp5ajWOr3sOPdAAXfCRfH0qnWAD8MlPO4c7zxvv+3t7SwWi+zt7b321+zu7qbv+3Pr4/HYawK4cnaggDuh67qMRqMza6PRKF3XDTQRcJfZgQLuhOc3ivsUHnAdBBRwZ0wmE8EEXAuX8AAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBo9EYBVZ76VCnl35dS/spVDwUAcJO96Q7Ux5N8JMk0yTullF99dSMBANxsbxpQ357kJ5J8a5KfTPJtVzYRAMAN9+ANv26c5LNJvjHJ+88ef1kp5Z0k7yTJw4cP8+TJkzefELh39vf380u/9Eu+dwA31psGVJ/kq5L81zzdjfq5F3+z1vrJJJ9Mkul0Wt96660PMSJwX/neAdxUb3oJ7yjJb6m1fibJdyT5D1c3EgDAzfamAfU4ydeVUn4qyU/XWn/hCmcCALjR3ugSXq21JvljVzwLAMCt4AdpAgA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQqNRa13uCUn4xyc+v9STAXfSxJJ8begjg3vmmWus3XHbQ2gMK4E2UUg5rrdOh5wC4iEt4AACNBBQAQCMBBdxUnxx6AICXcQ8UAEAjO1DArVVK+fqhZwDupwdDDwDcfaWU/5HkP39g+dcn+YO11n9bSnn/VV9fa91+9jzfnOTdWuufLKX8xiQ/mOSPXP3EAK8moIC1KqV8XZJ/k2QnyY8k+cFa638rpfzVJMtnh/3eJA+frf/JJJ+rtf6jUsrvSDJ/yVP/4SS/tZTy4y+s/cFa6/9ayx8E4AUu4QHr9qvyNIj+T5I/k+Svl1I+kmSU5FeeHfPRJP/ggq/94SRfSpJSyrcn+YE8jaYfTvK7knxLrfU7a63fmeSzz54TYO3sQAFrU0rZSfKnk3zsA5fp/mWSb07yHaWUv1Vr3S+l/LdSysdf+NqPJzmqtf7fJKm1HpVS+iTv5ul/3eBHk/zrUkqS/Lk8/X52uv4/FYCAAtboWRgt83SH6SDJoyR/s9Z6Wkr5e0l+qNb6c88O/+Ekn0/y/KeP/3KSv/6S5/27z375Q8/XSimjJF+8+j8FwHkCCli3X5Pk39Va/3sp5b8k+cellO9N8rVJfvnZLtXOC8c/zNN7o/5XkjzbYdqvte4/+/2vLaX8YK31r5VSviZP74XayNNI+8L6/zgAAgpYv1+X5HtKKZ9/9ngjT3ejvjXJP03y559/yi5JXryJ/MUnKaU8TPKn8nSH6h+VUn4kyW9K8qln//vuWqsdKOBaCChg3b45yW+rtX7pxcVSyk/XWr+j4Xm+IcnjJH+21roqpfxsrfVnnz1XSfKRK5sY4BICCli3P/piPD37BN13JvlMy5PUWo+THL/w+GdLKV+b5P1nS//kw48K8Hr8p1wAABr5OVAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQKP/B8hl7qNcVw82AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "rv = stats.norm()\n",
    "n_samples = 20\n",
    "ax.vlines(p_mean, 0, 21)\n",
    "for i in range(n_samples):\n",
    "    sample_ = samples[i]\n",
    "    s_mean_ = np.mean(sample_)\n",
    "    lcl = s_mean_ - rv.isf(0.025) * np.sqrt(p_var/n)\n",
    "    ucl = s_mean_ - rv.isf(0.975) * np.sqrt(p_var/n)\n",
    "    if lcl <= p_mean <= ucl:\n",
    "        ax.scatter(s_mean_, n_samples-i, color='gray')\n",
    "        ax.hlines(n_samples-i, lcl, ucl, color='gray')\n",
    "    else:\n",
    "        ax.scatter(s_mean_, n_samples-i, color='b')\n",
    "        ax.hlines(n_samples-i, lcl, ucl, color='b')\n",
    "ax.set_xticks([p_mean])\n",
    "ax.set_xticklabels(['母平均'])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:49.710540Z",
     "start_time": "2018-08-14T18:04:47.816297Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.951"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = stats.norm()\n",
    "cnt = 0\n",
    "\n",
    "for sample_ in samples:\n",
    "    s_mean_ = np.mean(sample_)\n",
    "    lcl = s_mean_ - rv.isf(0.025) * np.sqrt(p_var/n)\n",
    "    ucl = s_mean_ - rv.isf(0.975) * np.sqrt(p_var/n)\n",
    "    if lcl <= p_mean <= ucl:\n",
    "        cnt += 1\n",
    "cnt / len(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 正規分布の母分散(平均未知)の区間推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:49.714818Z",
     "start_time": "2018-08-14T18:04:49.711859Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 11.303,  15.767,  18.102, ...,  19.435,   9.265,  18.625])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_y = sample_u_vars * (n-1) / p_var\n",
    "sample_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:49.879900Z",
     "start_time": "2018-08-14T18:04:49.715894Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAloAAAFoCAYAAACPAVXRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlwleeB5/vvcxZtaEMSkkASCCR2EGBWm8UYb7GJk2AnTjudzk26PU6qu6bu3Kqpe9M9NXeqc2d6PDXVMzXVXXUdp9LXjl1px3ZC4o3YGLNYZt/RhkELAiSBEEgIgXS25/4h4RgswZE4R+9Zfp8qlescPeec3xOhk5+e9znva6y1iIiIiEjkuZwOICIiIpKoVLREREREokRFS0RERCRKVLREREREokRFS0RERCRKVLREREREokRFS0RERCRKVLREREREokRFS0RERCRKVLREREREosTjdICbCgoKbHl5eVRfw+fzkZKSEtXXiGXJPP9knjsk9/w19+ScOyT3/JN57jA+8z906NAla+2ku42LmaJVXl7OwYMHo/oaLS0tRLvMxbJknn8yzx2Se/6ae7nTMRyTzPNP5rnD+MzfGHMmnHE6dCgiIiISJSpaIiIiIlGioiUiIiISJSpaIiIiIlGioiUiIiISJSpaIiIiIlGioiUiIiISJSpaIiIiIlGioiUiIiISJSpaIiIiIlGioiUiIiISJSpaIiIiIlGioiUiIiISJR6nA4hI4in/6ftfua/lxY0OJBERcZZWtERERESiREVLREREJEpUtERERESiREVLREREJEpUtERERESiRJ86FJGw6JOEIiKjp6IlIo4Jt7yp5IlIvNKhQxEREZEoUdESERERiRIVLREREZEoUdESERERiRIVLREREZEoUdESERERiRIVLREREZEoUdESERERiRKdsFRExsVwJx0VEUl0WtESERERiRIVLREREZEouWvRMoN+YYw5bIz52QhjUo0xm40xR4wxPx6675+MMdVDX53GmJWRDi8iIiISy8LZo7UKSAOWAW3GmJestW23jXkaqAF+CDQZY16x1v5bGCxqwGfAwYilFhEREYkD4RStpcB2oAqoBpYAtxetpcAHDJayvUAlUDv0vfXAZ9baYATyikiEDbdJveXFjQ4kERFJPOHs0coBLgGFwI6h2yONsUD9bWP+N+BX95RSREREJA6Fs6LVA6QCdQyuXDWPMKYAOAY8M3QbY0wqsMBae2K4JzbGvAC8AFBSUkJLS8so449OV1dXVJ8/1iXz/JN57jD6+Yf7uxiN39lIv3Yy/+yTee6Q3PNP5rlDbM0/nKJ1CHjGWvuWMWYd8OoIY+YxuOJVBZweuv9hBvdnDcta+zLwMsCyZctseXl52MHHajxeI5Yl8/yTee5wp/nXfuWe4ceGO274seG699cO9zmTQzLPHZJ7/sk8d4id+Ydz6HAvkGuM2QccGGYjPMBmBkvVXuAVa+3A0P2bgC0RSSoiIiISZ+66omWttcDzN28bYx4EtgEV1tozQ2P6GSxVtz/230QuqoiIiEh8GcsleA4Ci/nqJw9FJEHocjkiIpEx6qJlre1j8JxZIpLkVMhERO5Ml+ARERERiRIVLREREZEoGcseLRERAAyWTDNAlhnAZ91ctyncwIvFOB1NRCQmqGiJJIDxuoxOjrlBhfsyua4b5Jh+sswAbmNvGROycMN66SOFjmAWjcF8um16xLOIiMQDFS0RuSODpczVzRzPRUrcvYSs4apNpcem0RrMpcem0WtT8RJkgvGTYXxMMD6yzAALPB1UeTu4HEqnMZhHczCfPpvi9JRERMaNipaIDMtDkDmeTua4L5Ll8tEX8nLIX8LJQAEDeMN6jjT8THdfYYani+Xe8yzznOdMaCIH/KVcs6lRnoGIiPNUtETkK/JNHw+mNJHjGqA9mMmBgTJaQ7mj3nvVj5f6YCH1wUKyTD8z3V3M81ygNLWbE4FiTgQmE9RnckQkgaloiSSor+7bqr3rvi2DZaGnnSWedm5YD38cmEV7KDsieXptGocDJZwMTmKZ5yxLvO1Uurs44C/jTCgXtIFeRBKQipaIAJBpBljnbabIfY2mwET2+Kfhi8JbRJ9NYae/gpPBq6z0nmVDaiOtwRw+9U2PyuuJiDhJ72oiQpmrm3UpTYBhp286TcE8or3C1BHK5p2Becx1X2C59zxPpdazzVcZ1dcUERlv2hwhkuSmuq6wIaWRHpvG7wfm0RTMZ7wO41kMdcFitvhm4TFBvp5aT319/bi8tojIeFDREklidXV1PJTSyKVQBh8OzKbPoU8CXgxl8W7/PLptOm+++Sbbt2/HWnv3B4qIxDgdOhRJIl/eID/dfZl13iY6Q5l85JtJALeDyeA6KWwZmM3/uD/Erl27uHjxIs888wwej96mRCR+aUVLJAnNcHexztvExRgpWTcFcfGNb3yDxx9/nIaGBt58800CgYDTsURExkx/Kookmenuy6z1NnMhlMXHvsqYKVk3GWNYtWoVXq+X9957jzfffBMXWYT0d6GIxCEVLZEkkmeus8bbzMVQJlt9lQRjrGR92dKlSwF477332JCSwye+CpUtEYk7KloiSSKVABtSTjOAh+2+ipguWTfdWrYax1S21r9UC9Tecl80LrgtIjIc/XkokgQMlgdTmsgwfj4ZqKA/zGsVxoKlS5fymW8aZe4eNqQ0Ygg5HUlEJGxa0RJJAks8bZS4r/KZbxqXbKbTce7oq5cOApgEPlidcob7va3s9k8b91wiImOhoiWS4Ka6rrDI287JQAGfByc5HWfMPg9OItM/wCJvBz02zek4IiJh0aFDkQSWbfpZm9JMZyiDff6pTse5Z4cDJTQHJrLcc46Ghgan44iI3JWKlkiCMoRYn9JICBfbfZUEE+LX3fCpfzqddgK/+93vaGtrczqQiMgd6dChSIJa6Okg33WDbQMV9NkUp+NETBAX2wYq+Tr1/K+XX+Hd/rlcJ3HmJyKJJRH+xBWR2+SYGyz2tNMcmEhraKLTcSKuHy8f+2biIcgjqafwEHQ6kojIsFS0RBKOZbW3hQAu9ibAvqyRdNt0dvgqyDM3WOVtdTqOiMiwVLREEsxcdydF7j72+cvi6nxZY3E+lMPRwGRmerqodF9yOo6IyFeoaIkkkAlmgKXec5wPZtMYzHc6zrg4FphCWzCL+72t5JobTscREbmFipZIwrA84D0DMHRCT+NsnHFiMezyzcCHi4dSGrVfS0RiioqWSIKocHdR6r7KIX8J12yq03HG1Q287PTNINv0c7/3DGCdjiQiAqhoiSQELwFWeM9xMTiBhmCh03Ec0RHK5mhgCpWey8zUfi0RiRE6j5ZIDBvuun8tL278yn1Vng7STICP/DOxSXLIcDjHA5MpcvWyyttKZyiTbpvudCQRSXJa0RKJcz09PczzXKAxkEeXneB0HEfd3K/lx826lCZchJyOJCJJ7q5Fywz6hTHmsDHmZyOMSTXGbDbGHDHG/HjoPo8x5g1jzAFjzH+KdHARGfTJJ58AcChQ4nCS2NCPl92+aeS7brDI0+50HBFJcuGsaK0C0oBlwAvGmCnDjHkaqAHWA/9gjEkFvgtstdYuB1YZY3IjE1lEbmpvb+f48ePUB4roS7IN8HfSGprIqUA+VZ52Csw1p+OISBILp2gtBbYDVUA1sOQOY1YBe4FK4FGgxBizA/iDtbY7EoFFZJC1lq1bt5Kens7xQLHTcWLOPn8Z120Ka1NacOsQoog4JJzN8DlAK1AI7Bi6PdyYS0AxUD90uxg4CmwAdhpj3rHWtn35QcaYF4AXAEpKSmhpaRnTJMLV1dUV1eePdck8/0Sa+83fk3PnztHc3MyKFSvw7VSRuJ0fD9X+cr6W+jlLvefY/6XLEUX7vSZWJNK/+7FI5vkn89whtuYfTtHqAVKBOgZXrppHGFMAHAOeGbp9DdhhrQ0ZY/YCM4Fbipa19mXgZYBly5bZ8vLysc1iFMbjNWJZMs8/Pude+5V71r9Ui8HyzdRa3KTyNzsD6HMtw2sPZVMXKGS+5yJng7m0h7KBeP23MDbJNNfhJPP8k3nuEDvzD+fd+RCw0lp7GlgHHBlhzDygk8FDjKcZPIS4cuj7N+8TkQiY6b7ERFc/B/2lhFSy7uigv4SeUCprvC14CTgdR0SSTDgrWnuBvzLG7AM+vP3w35DNwL8CPwB+aa0dMMb8HPi1MeaHwFvW2vORCi2SzNyEWOJt40JwAmdC+ozJ3QRxs8s3nY2pDSzznmePf9o9PV+45zYTEYEwipa11gLP37xtjHkQ2AZUWGvPDI3pBzbd9rhe4KmIphURZrk7yTB+dgRmkCzXM7xXl2wmdYEiFngv0BjMczqOiCSRsRxzOAgs5rb9ViISfS5CLPR00BHM5EIoy+k4ceVIYAq9oRRWe88QCOgQooiMj1EXLWttn7W2xlrrj0YgERlZpbuLCS4/xwKTnY4SdwK42eOfRq6rn+rqaqfjiEiS0C5akThhCFHlaaczlEHb0KfnZHTOh3JoDORRXV1NZ2en03FEJAmoaInEiRnuy2S5fBzzT0F7s8Zuv78Mr9fLe++9x+AWVBGR6FHREokDBssiTztdoXTOhoY7Z7CEqx8vjz32GK2trRw6dMjpOCKS4FS0ROJAufsKOa4Bjvkno9Wse7d48WLKy8v5+OOP6e3tdTqOiCQwFS2RmDe4mtUdSuNMaKLTYRKCMYavf/3rBINBPvzwQ6fjiEgCU9ESiXFTXd1MdN0Y+qShVrMiJT8/n9WrV1NbW0tTU5PTcUQkQaloicQ0yyJvO1dDqTTrRJsRt3r1aiZOnMiWLVsIBoNOxxGRBKSiJRLDJrt6KXBd53igGKvVrIjzer187Wtf49KlS+zdu9fpOCKSgMK51qGIOGSe5wI3rIemYL7TURLK7dcrfDglh507d7Jw4UKys3WOMhGJHK1oicSorq4uylw9nAxMIqhf1aja55+KtZaPPvrI6SgikmD07i0So/bv308IQ0NgktNREt41m8qaNWu0MV5EIk5FSyQG9ff3c+TIEZqDedwgxek4SUEb40UkGlS0RGLQkSNH8Pv91AUKnY6SNDwejzbGi0jEqWiJxJhQKMS+ffuYOnUqXXaC03GSyqxZs5g1axa7du3i2rVrTscRkQSgoiUSY06ePElPTw+rVq1yOkpSeuyxxwgEAnzyySdORxGRBKCiJRJj9u7dS25uLrNnz3Y6SlLKz89n5cqVHDlyhPb2dqfjiEicU9ESiSFtbW20trayYsUKXC79ejpl3bp1ZGRk8Mc//hFrrdNxRCSO6Z1cJIbs27ePlJQUlixZ4nSUpJaWlsaGDRtobW2lrq7O6TgiEsdUtERixLVr16ipqWHRokWkpaU5HSfpLVmyhKKiIrZu3Yrf73c6jojEKV2CRyRGHD58mFAoxIoVK5yOkpRuvywPwI6ffI1XX32VPXv2sG7dOgdSiUi804qWSAyw1nL48GHKy8spKChwOo4MKS8vZ+7cuVRXV3P16lWn44hIHNKKlogDbl89KXH18FhqD4888ohDiWQ45T99n0zjYVNqgL/577/iU//0EccNp+XFjdGMJyJxQCtaIjFgtqeTjIwM5s6d63QUuc01m0pdoIhKTxf5ps/pOCISZ1S0RByWjo8yVzeLFy/G7XY7HUeGcTxQTL/1sNx7DtDpHkQkfCpaIg6b5bmEy8DSpUudjiIj8OPhiH8Kk929lLl6nI4jInFERUvEQQbLLPcl2oJZ5OXlOR1H7uBksIDuUBrLvWcxhJyOIyJxQkVLxEElrh4yXT5OBiY5HUXuwuLioL+UHNcAs92XnI4jInFCRUvEQbM9l7hhPbSGcp2OImE4G8qhPZjFEm8bXgJOxxGROKCiJeKQDHyUuro5FSggpF/FOGE44C8llQBVng6nw4hIHNB5tEQcMnNoE/znwcETlI50LiaJLV12Ao3BfOZ5LnAyOIlrNtXpSCISw/RntIgDbm6CPx/Mptfquobx5nCgBID7POcdTiIise6uRcsM+oUx5rAx5mcjjEk1xmw2xhwxxvx46L6/McacMMZUG2O2RDq4SDzTJvj41mdTqA0UU+G5rJOYisgdhbOitQpIA5YBLxhjpgwz5mmgBlgP/IMxJhXIBP6ttXaNtfaJCOUVSQiVni76rYezoRyno8gYndBJTEUkDOEUraXAdqAKqAaW3GHMKmAvUMlg0fo7Y8x+Y8zfRyauSPy7fv06U13dNAbztAk+jvlxc9Q/mcnuXkpcuuC0iAwvnM3wOUArUAjsGLo93JhLQDFQP3S7HvjIWvupMWa7MWaatfbMlx9kjHkBeAGgpKSElpaWMU4jPF1dXVF9/liXzPOPpbnX19fjNpZTgQKno8g9OhmcxLzQRZZ7z9E2kI3F3PL9aL+n3U0s/bt3QjLPP5nnDrE1/3CKVg+QCtQxuHLVPMKYAuAY8AzQY6399Ze+fwCYDtxStKy1LwMvAyxbtsyWl5ePMv7ojcdrxLJknn+szP2jjz6iK5TOFZvhdBS5RyFcHPKX8FBqExXuLk4Hby3PsfBvLhYyOCmZ55/Mc4fYmX84xy0OASuttaeBdcCREcbMAzoZPMR42hjzz8aYhcYYF4P7u05FKLNI3Lpw4QLt7e1azUogLaGJdIYmcJ/3PG5dmkdEbhNO0doL5Bpj9gEHrLVtw4zZDDw8NPYVa+0A8E/A/8tgCfvQWqvPQUvSO3bsGC6Xi6agrmuYOAZPYjrB+JnnueB0GBGJMXc9dGittcDzN28bYx4EtgEVN/dcWWv7gU23Pe4ksCaiaUXiWDAY5Pjx48yaNYuBI16n40gEXQhl0RrMocrTweeBAgbQz1dEBo3lI08HgcXAcCtbIjKCxsZG+vr6WLx4sdNRJAoO+kvxEGSRt93pKCISQ0ZdtKy1fdbaGmutPxqBRBLV0aNHycjIoLKy0ukoEgU9Np1TwQLmuDvJNANOxxGRGKGT+IiMg+vXr3Py5Emqqqpwu91Ox5EoOeqfgsXo0jwi8gVdVFokysp/+j5z3RdZlRLipzuu8uPtunh0orpOCrWBQhZ5O6gJFA97ofCWFzc6kExEnKIVLZFxUOm5xKVQhs6dlQRqAsX0WzfLvOecjiIiMUBFSyTKJprrFLiuczqQ73QUGQc+PBwPTKbEfZXJujSPSNJT0RKJsgp3F0FrdO6sJNIQKORaKGVoVUsXnBZJZipaIlEUCoWY4bnM+VCOzq2URIK4OByYQoHrOtPdV5yOIyIOUtESiaIzZ84wwfhp1GpW0mkK5nM5lM59nvO4dGkekaSloiUSRcePH8dnXZwN5jodRcaZxXDIX0K2a4BZ7ktOxxERh6hoiUSJ3++nvr6eM8GJBPWrlpTOhXJoD2ay2NuGh6DTcUTEAXr3F4mSzz//nIGBARqD+rRh8jIc9JeSbgIs0AWnRZKSipZIlJw4cYLMzEw6QllORxEHXbKZtARzWeDpIA1duUwk2ahoiUTBjRs3OHXqFAsWLMBinI4jDjvsL8FNSBecFklCKloiUVBbW0soFKKqqsrpKBIDbl5wera7kytXdLoHkWSiax2KRMGJEycoKCiguLjY6SgSI476p1Dh7uJ//x+vs8s/45bv6fqHIolLK1oiEdbd3U1raysLFy7EGB02lEHXSaEuUESF5zJ55rrTcURknKhoiUTYiRMnAFi4cKHDSSTWnAgUM2DdLNUFp0WShoqWSARZazlx4gRlZWVMnDjR6TgSY25ecLrUfZViXXBaJCmoaIlE0IULF+js7NRqloyoPlBIX8jLMu95dMFpkcSnoiUSQcePH8flcjF//nyno0iMGrzgdAmTXH1Mc+kTiCKJTkVLJEKstdTW1lJRUUFGRobTcSSGNQbzuRJKY6n3PEYXnBZJaCpaIhFy9uxZrl69yoIFC5yOIjFu8ILTpeS4Bpjl7nI6johEkYqWSITU1NTg8XiYPXu201EkDpwN5dAxdMFpn8/ndBwRiRKdsFQkAkKhEHV1dcyaNYvU1FSn40hcMBwKlLIxtYG9e/eybt26uz6i/Kfvf+U+nexUJLZpRUskApqbm+nr69NhQxmVi6FMzgRz+eyzz7h+XScxFUlEKloiEVBTU0NKSgqVlZVOR5E4c8hfgt/vZ9euXU5HEZEoUNESuUeBQICGhgbmzp2L1+t1Oo7EmR6bzuLFizlw4IAuOC2SgLRHS+QeNTY20t/fz/z584fdQyNyN+vXr+fEiRPs2LGDTZs2AcPvxxKR+KMVLZF7VFNTQ3p6OjNmzHA6isSp7OxsVq5cyfHjx+no6HA6johEkIqWyD3w+XycPHmSefPm4Xa7nY4jcWzNmjWkp6fz8ccfOx1FRCJIRUvkHnz++ef4/X592lDuWVpaGmvXrqWxsZGmpian44hIhKhoidyD2tpasrKymDp1qtNRJAEsX76cnJwctm7dii44LZIY7lq0zKBfGGMOG2N+NsKYVGPMZmPMEWPMj2/73j8bY96IVGCRWNHf38+pU6eYN28eLpf+ZpF75/F42LBhAx0dHcxwX3Y6johEQDj/77AKSAOWAS8YY6YMM+ZpoAZYD/yDMSYVwBjzKJAdmagisaWhoYFgMKjDhhJRCxcupLi4mPs853HpgtMicS+corUU2A5UAdXAkjuMWQXsBSqNMbnA3wD/T2SiisSW2tpacnNzKSkpcTqKJBBjDI888ghZLh9z3BedjiMi9yic82jlAK1AIbBj6PZwYy4BxUD90O1/B/wU8I/0xMaYF4AXAEpKSmhpaQk/+Rh0dXVF9fljXTLPP9Jz7+/vp7Gxkfnz53PmzJmIPrckn9vf+9xuN+eD2SzytnM6WIDvDm/Vd3rfTObfeUju+Sfz3CG25h9O0eoBUoE6BleumkcYUwAcA54B+oA1wEwGDzvOMMb8hbX2tS8/yFr7MvAywLJly2x5efnYZjEK4/EasSyZ5x/JuR86dAhrLQ888ACTJ0/+0ndqI/YakjyG+7d50F/KN1LrWOjp4FCgdFSPHc33E10yzz+Z5w6xM/9witYh4Blr7VvGmHXAqyOMmcfgilcV0GCtnQtgjCkHXry9ZInEs7q6OvLy8iguLnY6iiSA4c8Cn0FjMI95ngs0BAvpsynjnktE7l04e7T2ArnGmH3AAWtt2zBjNgMPD419xVo7EMGMIjGlr6+P5uZm5s+fjzHG6TiSwA4HSjDAEs95p6OIyBjddUXLWmuB52/eNsY8CGwDKqy1Z4bG9AObRnh8C/BnkQgrEgvq6+ux1jJ//nyno0iC67Op1AUKWeC5QF2giMs2w+lIIjJKYzn5z0FgMTDcypZIwqutraWgoIDCwkKno0gSOB6YzABulnnPopOYisSfURcta22ftbbGWjvipwlFEtW1a9doaWnRYUMZNz48HPNPocTdS4nrqtNxRGSUdDprkVGoq6sD0GFDGVcNwUlcDaWy3HsOo1UtkbiioiUyCrW1tRQWFjJp0iSno0gSCeHioL+Eia4bVLovOR1HREZBRUskTFevXqW1tZV58+Y5HUWS0JnQRC4GJ3Cftw0PQafjiEiYVLREwqTDhuIsw4FAGRnGz3zPBafDiEiYVLREwlRbW0tRUREFBQVOR5EkdTGUSUtwIgs9HaSPfHUzEYkhKloiYejp6eHcuXNazRLHHfSX4MayxKuTmIrEAxUtkTDU1g5ew1BFS5zWa9OoD05ipvsSE811p+OIyF2oaImEoba2lsmTJ5OXl+d0FBGO+qfgx81y71kGL94hIrFKRUvkLq5cuUJbW5s+bSgxw4eHI0MnMT116pTTcUTkDlS0RO5CnzaUWNQQnERPKJWPPvqIYFCnexCJVSpaIndRW1vLlClTmDhxotNRRL5gcXHAX0ZXVxcHDx50Oo6IjMDjdACRWHb58mXa29t59NFHv/K98p++70AikT85G8ph+vTp7Ny5k6qqKtLT052OJCK3UdESuYObhw21P0tik+Gxxx7j5z//Obt27eLxxx8fceRwfxi0vLgxmuFEBB06FLmj2tpaSktLyc3NdTqKyLCKi4tZsmQJ+/fvp6ury+k4InIbFS2REXR1ddHR0aHVLIl5GzZswOPxsHXrVqejiMhtdOhQZAQ6bCjxIjMzk7Vr17Jt2zaamprY8HL9l75b61guEVHREhlRbW0tZWVlLPqv1U5HEbmrVatWcfjwYf74xz9imIrFOB1JRNChQ5FhXbp0iQsXLmg1S+KGx+Phscceo7OzkznuTqfjiMgQFS2RYdy8tqGKlsST2bNnM2PGDJZ4z5NKwOk4IoKKlsiw6urqmDp1KtnZ2U5HEQmbMYbHH38cL0GWeM87HUdEUNES+YrOzk4uXryo1SyJS4WFhTQEC5nt7mSiue50HJGkp6IlchsdNpR4d8Q/BR9uVnrPAtbpOCJJTUVL5DZ1dXVMmzaNrKwsp6OIjIkPD0f8JUx29zLN1e10HJGkpqIl8iUXL16ks7NTq1kS904GJ3EllM5y71nchJyOI5K0VLREvqSmpgZjjIqWxD2LYa+/jCyXj4WedqfjiCQtFS2RIdZaamtrKS8vJzMz0+k4IvesI5RNUyCPhZ4Osky/03FEkpKKlsiQjo4OLl++zPz5852OIhIxB/ylWAwrvGedjiKSlFS0RIbU1tZijGHu3LlORxGJmOukcDQwhanuHsq0MV5k3KloifCnw4YzZswgIyPD6TgiEVUbKKQ7lMZKb6s2xouMMxUtEaCtrY3u7m4dNpSEZHGx1z91aGN8h9NxRJKKipYIg4cNXS4Xc+bMcTqKSFS0h7JpDkxkoaedTDPgdByRpHHXomUG/cIYc9gY87MRxqQaYzYbY44YY348dN8iY8xuY0y1MeZVY4yJdHiRSLh52LCyspL09HSn44hEzX5/GRYzdMZ4ERkP4axorQLSgGXAC8aYKcOMeRqoAdYD/2CMSQXcwNPW2jXAZKAoIolFIuzcuXNcvXpVhw0l4Q1ujJ/MVHe3NsaLjBNPGGOWAtuBKqAaWAK0DTPmAwZL2V6g0lp72BjzXWPMfwZ+a63VxgCJSTU1NbjdbmbPnk35T993Oo5IVNUGiqh0d7HK24rP5yMlJcXpSCIJLZyilQO0AoXAjqHbw42zj6y8AAAgAElEQVS5BBQD9TfHWGt/Y4zZDLxpjJlmrT0TidAikRIKhairq2PmzJmkpqY6HUck6iwudvunsTH1JN/72b9wMFB2y/dbXtzoUDKRxBRO0eoBUoE6BleumkcYUwAcA54Beowx9wHHrbU+Y8xOBle7bilaxpgXgBcASkpKaGlpGeM0wtPV1RXV5491yTz/kebe0dHBtWvXKCwsjPq/P5FoGMu/24uhLE4GCpjvuUBjMJ8r9k+nNEmk3wO95yWvWJp/OEXrEPCMtfYtY8w64NURxsxjcMWrCjgN/C/gNwwedlwCfHb7g6y1LwMvAyxbtsyWl5ePfgajNB6vEcuSef7Dzb22thaPx8Pq1auHDqHUjnsukXsx/O/03f8dH/KXMs3dzQMpZ3h/YA5g7vB88SvR5jMayTx3iJ35h7MZfi+Qa4zZBxyw1t6+PwtgM/Dw0NhXrLUDwH8G/s4Y8xnQbq3dH6nQIpEQCoWor69n1qxZ2qciSWcAD/v9pRS6+pjlvuR0HJGEddcVLWutBZ6/edsY8yCwDai4uefKWtsPbLrtceeARyOaViSCmpub6evr06cNJWk1BvOZGeximfccrcFc+vE6HUkk4YzlhKUHgcV89ZOHInGlpqaGlJQUZs6c6XQUEYcY9vin4iHEcp1bSyQqRl20rLV91toaa60/GoFExkMgEKC+vp65c+fi9eqveElePTadE4FiKj2XmeLqcTqOSMLRJXgkKZ0+fZqBgQEWLFjgdBQRxx0PTKYnlMoD3jP4fD6n44gkFBUtSUo1NTVkZGQwffp0p6OIOC6Ii8/85WS5fGzfvt3pOCIJJZzTO4gkFJ/Px8mTJ1m8eDFut9vpOCL3JFJXM7gQyqIhMAmzbx8LFiygpKQkIs8rkuy0oiVJp6GhgUAgoMOGIrc56C8lMzOTd955h2Aw6HQckYSgoiVJp6amhuzsbKZOnep0FJGY4sfNxo0buXjxIp999pVzTIvIGKhoSVK5fv06jY2NzJ8/H2OM03FEYs7s2bOZP38+u3btorOz0+k4InFPe7QkqdTX1xMKhVi4cKHTUURi1hNPPEFTUxPvvvsuP/rRj5j+tx98ZYwuPi0SHq1oSVKpqakhPz+f4uJip6OIxKwJEybw+OOPc/bsWfbv19XTRO6FipYkjatXr9LS0sKCBQt02FDkLqqqqpg5cyYff/wx2abf6TgicUtFS5JGbW0tgD5tKBIGYwxPPfUUHo+HNd4WDNbpSCJxSUVLkkZNTQ3FxcUUFBQ4HUUkLmRlZfHEE09Q5L7GPM8Fp+OIxCUVLUkKV69epa2tTatZIqO0cOFCzgRzuc9znhxzw+k4InFHRUuSQlNTE6DDhiKjZYxhj28aAVysSdEhRJHR0ukdJOFZa2lqamL69Onk5OQ4HUckpg1/SR8ve/3TWJ/SxAJPBycCk8c9l0i80oqWJLxz587R29urc2eJ3IPm4ERaghNZ4mkjV4cQRcKmoiUJ7/jx47jdbubNm+d0FJE4Ztjjm4oPN+tSmggEAk4HEokLKlqS0ILBILW1tZSVlZGamup0HJG41o+Xal85+a4bbN++3ek4InFBRUsS2qlTp7hx4wYVFRVORxFJCOdCuTQEJrF7926am5udjiMS87QZXhLaiRMnyMjIYMqUKV/53vCbfkXkbg74S1ldFOL3v/89P/nJT0hPT3c6kkjMUtGShNXf38/JkydZunQpG16uB+qdjiSSEAK4efrpp/nlL3/J+++/zzPPPKPLWomMQIcOJWHV1dURDAapqqpyOopIwpkyZQrr16+ntraWEydOOB1HJGapaEnCOn78OPn5+cMeNhSRe7d69WrKysr44IMP6O7udjqOSExS0ZKE1N3dzZkzZ6iqqtIhDZEocblcbNq0CWstv/3tbwkGg05HEok52qMlCenmoQwdNhSJji9/mGS6u5T1via2b9/OI4884mAqkdijFS1JONZajh8/ztSpU8nNzXU6jkjCaw7mcTJQwGeffcbp06edjiMSU1S0JOG0t7dz6dIlrWaJjKN9/jIKCwvZvHkzvb29TscRiRkqWpJwjh49qkvuiIyzIG6+/e1v4/f7+d3vfkcoFHI6kkhMUNGShBIIBDhx4gRz587VSRRFxtmkSZN44oknaGlp4dNPP3U6jkhMUNGShNLQ0EB/fz9LlixxOopIUlq8eDFVVVXs3LmTlpYWp+OIOE5FSxLK0aNHycnJYfr06U5HEUlKxhiefPJJ8vLy+O1vf6v9WpL0VLQkYfT09NDY2MjixYt17iwRB6WmpvLss88yMDDA22+/rfNrSVK7a9Eyg35hjDlsjPnZCGNSjTGbjTFHjDE/HrpvujFmmzHmoDHm7yMdXOR2R48eBQYPXYiIswoLC3nqqadobW1l69atTscRcUw4JyxdBaQBy4A2Y8xL1tq228Y8DdQAPwSajDGvAP8B+Ftr7f6hsvWP1tqrkYsu8ifWWo4ePcr06dN17iwRh3z5JKY3vbR+Jfv27aO0tJQFCxY4kErEWeEcOlwKbAeqgGpguF3GN8esAvYClcBvgMND3+8GBu41rMhIWlpa6O7u1iZ4kRjz6KOPUlZWxjvvvMPFixedjiMy7sIpWjnAJaAQ2DF0e6QxFqgHcqy1W4GgMea/AL+21qpoSdQcPXqU1NRU5syZ43QUEfkSt9vNd77zHVJSUnjzzTfp7+93OpLIuArn0GEPkArUMbhy1TzCmALgGPDM0G2AfwJ2WWvfHO6JjTEvAC8AlJSURP2jwF1dXVF9/liXqPP3+XzU1tZSUVHB+fPnnY4jIl9y8319zZo1fPTRR/z617/moYceGpcPrCTqe144knnuEFvzD6doHQKesda+ZYxZB7w6wph5DK54VQGnjTHfBc6OVLIArLUvAy8DLFu2zJaXl48u/RiMx2vEskSc/6FDhwgGg6xbt46SkpIRRtWOayYRGXTzPae8vBxjDH/84x9pbm5mw4YN4/r6ySiZ5w6xM/9wDh3uBXKNMfuAA8NshAfYDDw8NPaVocOEfwM8a4ypHvpaGLHUIl9y5MgRCgsLmTJlitNRROQOVqxYwZIlS/j00085ceKE03FExsVdV7SstRZ4/uZtY8yDwDagwlp7ZmhMP7Dptseti2xUka9a/LdvsyntPPt8pfy3v/0AgJYXNzqcSkSGY4xh48aNdHV18c4775CXl3eHVWiRxDCWE5YeBBYDw61siYyrWZ5OgtbQFMx3OoqIhMHtdvPss88yYcIEfvOb3+jM8ZLwwtmjdQtrbR+D58wScZTf76fS3cWZYC79eL+4f7hz+YhI7JgwYQLPPfccv/zlL3njjTf44Q9/iNfrvfsDReKQLsEjcaumpoZUE6QhWOh0FBEZpaKiIp5++mna2tp45513GNylIpJ4VLQkbh08eJAroTQuhDKdjiIiYzBnzhwefvhhampq+OSTT5yOIxIVoz50KBIL2traaGtr42RgKqALSIvEq9WrV3PlyhWqq6vJzc1l6dKlTkcSiSgVLYlLBw4cwOv1cvpGntNRROQe3PwkYm9vL++//z7Z2dnMnDnT6VgiEaOiJXHnxo0b1NTUUFVVhf8zHf0WiXcul4tvf/vbvPLKK7z11lv86Ec/YvLkyaN6juE+BLPjJ/MjFVFkzPT/UhJ3jh07RiAQYNmyZU5HEZEISUlJ4bnnniMjI4Nf//rXdHd3Ox1JJCJUtCSuWGs5dOgQJSUlo/6LV0Ri28L/sov/70IpXb03+Pv/+XNm//T3TkcSuWcqWhJXWlpauHTpklazRBJUj01nm6+STDPAo6mnGBgYcDqSyD1R0ZK4cvDgQdLS0pg/X3svRBLVhVAWO3wV5JvrvPHGGwQCAacjiYyZipbEjd7eXhoaGliyZInOIi2S4M6GcvnUP52WlhbefvttQqGQ05FExkRFS+LG4cOHCYVCOs+OSJJoCubzxBNPcPLkSZ09XuKWTu8gcSEQCHDgwAEqKyvJz9cFpEWSxYoVK+jv72f79u2kpqbyta99DWN0kmKJHypaEhdOnDhBX18f999/v9NRRGScrV27lhs3brB37148Hg//ZtsAuiKExAsVLYl51lr27NlDUVER06dPdzqOiIwzYwyPPfYYgUCA3bt3s8xTzMFACSpbEg9UtCTmNTY20tnZybe+9S0dMhCJM8OdsR2g5cWNo3oeYwxPPvnk4I2DB7HAIZUtiQPaDC8xb8+ePWRlZbFgwQKno4iIg26WrYbAJKq8HdznOQ9og7zENq1oSUzr6OigqamJhx9+GLfb7XQcEXGYMYY9/qkYLIu8HYDhcGAKWtmSWKWiJTFt7969eL1endJBRL7EsNs/DYBF3nZchDgYKEVlS2KRipbErN7eXk6cOMHSpUtJT093Oo6IxJTBshXCsNB7gRQTZI9/GlZlS2KMipbErP379xMKhVi1apXTUUQkJhn2+qfis24WeTvwEmKXvxyr7ccSQ1S0JCb5fD4OHjzI3LlzycvLczqOiMQsw+FAKT7cLPeex2uCbPdVEFTZkhihoiUx6ejRo/T3999ygtKRPiYuIlITmIzfurnf28qjKZ+zzTfT6UgigIqWxKBgMMju3bspLS2lrKzM6TgiEgXR+MPpZLAQP27Wept5IrWBvr6KiL+GyGhpbVViztGjR+np6eHBBx90OoqIxJmmYD5bfTPJNgNs2bKFixcvOh1JkpyKlsSUYDDIp59+SklJCRUV+mtUREavLZTDBwNzCIVC/Mu//AvNzc1OR5IkpqIlMeXmatb69et1uR0RGbPLNoMnn3yS7OxsXn/9dY4fP+50JElSKloSM7SaJSKRlJmZyY9+9CPKysrYvHkzO3fuxFpdskfGl4qWxAytZolIpKWnp/P973+fqqoqduzYwdtvv43P53M6liQRfepQYoJWs0RkOJH4dKLH4+Fb3/oWRUVFfPzxx3R1dfFnf/Zn5ObmRiChyJ1pRUtiglazRCSajDE88MADfO9736O7u5tf/OIXnDlzxulYkgS0oiWO02qWiIyXyspKnn/+ed544w1+9atf8fjjj/Ps7y5y+wWpW17c6ExASTh3XdEyg35hjDlsjPnZCGNSjTGbjTFHjDE/HrrPbYz5B2NMdaRDS2LRapaIjKeCggKef/55Kioq2LJlCw96m/EQdDqWJKhwVrRWAWnAMqDNGPOStbbttjFPAzXAD4EmY8wrgAU+AzZELK0kHL/fz65du76ymqXL7YhINKWlpfHcc89RXV1NaNsn5Lv62O6r4IrNcDqaJJhwitZSYDtQBVQDS4Dbi9ZS4AMGS9leoNJaWwu8b4z5j5GLK4lmz549XL16laefflqrWSIyrowxrF27ln+/5TzrU5r4emoDe/xTOR0sGPaPPR1OlLEIp2jlAK1AIbBj6PZwYy4BxUD9CGNEbtHb20t1dTVz585l2rRpTscRkQSz/qVaoPaW+4YrSxdCWfyhfx4PpjSxNqWFyYFe9vqn4sc9TkklkYVTtHqAVKCOwZWr4a5l0AMUAMeAZ4Zu35Ux5gXgBYCSkhJaWlrCediYdXV1RfX5Y12szb+6uppgMMicOXOi/rMXEQFGfK/px8tHvlks8rSxyNNOkauXXf7pXAxl3fWxsSjW3u/HWyzNP5yidQh4xlr7ljFmHfDqCGPmMbjiVQWcDufFrbUvAy8DLFu2zJaXl4fzsHsyHq8Ry2Jl/m1tbTQ2NvLAAw9QVVU1zIjaYe4TEbk3w78HDr7fWAxHAyW0hXJY523iiZSTHA9M5mhgMhZXzLx/hive8kZarMw/nPNo7QVyjTH7gAPDbIQH2Aw8PDT2FWvtQAQzSoKx1vLhhx+SkZHB2rVrnY4jInKLi6FM/jAwn9PBfBZ729mY2kC26Xc6lsSpu65o2cELQz1/87Yx5kFgG1BhrT0zNKYf2DTC41dFJqokivr6elpbW9m4cSNpaWlOxxER+Qo/bj7zT+dcMIfVKWf4Zmotu3dPZtWqVbhcOte3hG8s/1oOAov56icPRe4qEAiwdetWCgsLue+++5yOIyJyR2dCeWzun09bKIetW7fyy1/+kgsXLjgdS+LIqM8Mb63tY/CcWSKjtm/fPrq7u/n+97+vvwpFZNyN5Rx9N0hhm6+C975fzpYtW3j55ZdZu3Yta9euxe3WJxPlznQJHhk3PT097Nq1i1mzZulSOyISZwxff/0MqcxkhfcsoZ07+f32ffy7H303ZjZdS2zSkoKMC2st7777LtZavva1rzkdR0RkTAbw8ql/BlsHKvEQ4tVXX+V3v/sdvb29TkeTGKWiJePi2LFjNDY28sgjjzBx4kSn44iI3JNzoVw2D8xn3bp11NXV8c///M/s3buXUCjkdDSJMSpaEnW9vb18+OGHTJ06leXLlzsdR0QkIoK4eeihh/jrv/5rpk6dyocffsjPf/5zGhsbnY4mMURFS6LKWst7771HIBDgG9/4hq5nKCIJJy8vj+9973t897vfxe/38/rrr/P666/r04kCaDO8RFlNTQ2ff/45jz32GPn5+U7HERGJCmMMc+bMobKykgMHDrBr1y5eeuklFi9ezEMPPUR2drbTEcUhKloSNdeuXWPLli2UlpaycuVKp+OIiESdx+Ph/vvvZ/HixXz66afs37+fmpoali9fzurVq5kwYYLTEWWcqWhJVFhr+eCDD/D5fHzzm9/84pxZw53DpuXFjeMdT0QkqtLT03nsscdYvnw5O3bsYM+ePRw8eJCVK1fywAMPkJ6e7nREGScqWhIVhw4dor6+nocffpiCgoI7jh3LCQRFROLBxIkT2bRpE2vWrGHnzp1UV1ezf/9+Vq1axapVq1S4koCKlkTcuXPn2LJlCzNnzmT16tVOxxERcdykSZP49re/zdq1a9m5cye7du1iz549LF26lPvvv197uBKYipZEVF9fH2+++SbZ2dls2rRJnzIUEfmSoqIinn32WS5cuMDu3bvZt28f+/fvp6qqigceeIBJkyY5HVEiTEVLIiYUCvH2229z48YN/uqv/kpL4iIiIygqKmLTpk089NBD7N69myNHjnD06FFmzpzJypUrmTFjhv5QTRAqWhIx27Zto6WlhW9961sUFxc7HUdEJObl5uby5JNP8uCDD7J//34OHTrE66+/TkFBAStWrGDRokWkpKQ4HVPugYqWRERdXR27d+9m2bJlLFq0yOk4IiKOCffT1cONczGbzc9NZf/+/XzwwQds27aNRYsWcd9991FUVBSVvBJdKlpyz9rb2/nDH/5AaWmpLhgtIkkl3E9NhzsuhItFixZRVVXF+fPnv1jl2r9/PyUlJSxZsoQFCxaQmpp6L7FlHKloyT25ePEir732Gunp6XznO9/B7XY7HUlEJO4ZYygtLf3iD9jjx49z+PBh3nvvPT788EPmz5/PwoULKS8v/+I8hRKbVLRkzC5fvsxrr72G2+3mBz/4gT6eLCISBRkZGaxatYqVK1dy/vx5Dh8+TF1dHUePHiUrK4sFCxawcOFCiouLtYE+BqloyZh0d3fzq1/9ilAoxA9/+EPy8vKcjiQiktC+vMr1xBNP8Oh/eoOK4GV6du9lz549dIfSeOrB5cybNw9rrdNxZYiKloxab28vr732GgMDA/zgBz/QeV9ERMaZ1+ulJZRHiy+PVAKUuy8z3X2F6upqPv30U7KysqiqqmLu3LlMmTJFK10OUtGSUbl27RqvvfYa165d4y/+4i+YPHmy05FERBLKaK8JO4CHk8FCTgYLqf2/1tPQ0MChQ4fYvXs3n332Gdetl3PBHM4Gc2gLZRPArWvMjiMVLQlbR0cH//qv/8r169f58z//c0pLS+84XtcwFBGJjHDfTydMmMDSpUvJz8+nsLCQU6dO8T/f3km5+wqzPJcIWENHKIs9e/KYMWMGhYWFWu2KMhUtCUt9fT2bN28mPT2dv/zLv9RKlohIDLq1kNUO/bcC4w9R5LpGmauHUncPH330EQCZmZnMmDGDGTNm8L03Gumzt542Qitf905FS+7IWkt1dTWffPIJJSUlfPe73yUrK8vpWCIiMgoWFx2hbDpC2RwIlHHsb9fQ1NREU1MTp0+f5vjx4zybBr2hFDpCWXSEsrgQysRaqxWve6SiJSPy+/28++67nDhxgoULF/LUU0/h9XqHHavDhCIi8SMnJ4clS5awZMkSrLVcuHCBH/3TBxS7eylz9zDT0wXAP/5jM2VlZV982nHKlCkj/v+ADE9FS4bV0tLCe++9R1dXFxs2bGDNmjX6q0ZEJAEZYyguLqY+WER9sAiw5Jp+ily9/KQql3PnztHQ0ACAy+WisLCQyZMnM2XKFKZMmUJhYSEej+rESPS/jNyiv7+frVu3cvjwYXJzc/n+979PRUWF07FERGTcGLptOt3BdDZtGtyj1dfXx7lz5zh37hxtbW3U19dz5MgR4E/lq6ioiKKiIoqLiykqKiIjI8PJScQMFa0kNNJHhxsaGnj//ffp6+vj/vvvZ/369bpqvIhIAhrtdo8JEyYwe/ZsZs+ePfTYiWQaHwWmj3zXdfLOXyevvY4Mc+yLx2RmZjJp0qRbvgoKCsjIyEiqIyQqWknPUuQaPDdWU1MTRUVFPPfcc0yZMgUY/flcREQkcYxcyAzXbCrXbCotoT9dGSQNPxNdN8hzXeeFygI6Ozs5evQoPp/vT2PS0sjPzyc/P5+8vDzy8/OZOHEiEydOJD09PeFKmIpW0rKUunqo8rRT5O7jwoUJPProo6xcuVIXhhYRkTHpx0t7yEt7KJtvfnPwj3JrLVevXqWzs5NLly7R1dXF5cuXOXPmDMePH7/l8ampqV+UrpycHHJycsjNzf3iv2lpaXFXxFS0ksyNGzeodF9inucC+a4b9IZS2OObyrt/9319kkRERCLOGPNFaaqsrLzle36/n8uXL9Pd3c3ly5e5cuUK3d3ddHZ2curUKQKBwC3jPR4P2dnZX3xlZWWRlZVFZmbmF//NzMwcz+ndlYpWEhgYGODkyZMcOHCAtrY21qaE6A6lsctXTlMwD4uLmf/xo7CfT6dyEBGRSPB6vV9sor+dtZbr16/T09NDd3c3PT09XL16ld7eXq5evcqZM2fo7e0lFAp95bH33Xcf5eXl4zCDu7tr0TKDa3QvA0uB96y1//cwY1KBN4By4CVr7c+NMXnAW0A+8B+tte9GMriMzO/3c/78ec6ePUtrayvNzc0Eg0EyMjJYsWIF/2FHN5fsBCC+ll9FRCR+3OseX2MMEyZMYMKECV/sG76dtZY5f/sHMoyfdOP/4r9PDFPcnBLOitYqIA1YBrQZY16y1rbdNuZpoAb4IdBkjHkF+DHwa2Af8DagohVhoVCInp4eurq6vvg6f/48HR0dXzT8goICli5dyoIFCwgEAkyfPp0fb9eKlIiIjL9wj4iEW8iMMQzgZcB6uWL/dH9hYeFY4kVFOEVrKbAdqAKqgSXA7UVrKfABg6VsL1A5dN//yWBBu2yMmWCt7YtQ7rhlrcVaSygUwlpLMBi85SsQCODz+W75GhgYoK+v75ava9eu0d3dTTAY/OK5U1JSmDx5Mg888ABlZWWUlZUx9+8/gbMh2Hlzw2GdMxMXEREJUyJtUQmnaOUArUAhsGPo9nBjLgHFQP3Q7Zv3XWGwmOUAjhWtzz//nLfeeguXyzWmx1tr7z7otrFf/u/Nr3uRlpbGhAkTyMzM5ODFENdCk+ixqVy1aVwNpXHjhgd6DDT0A6eGvkRERMQp4RStHiCVwaWQpUDzCGMKgGPAM0O3e4AKBlfB/v3Q7VsYY14AXhi6ec0Yc3KU+UergMHyl6ySef7JPHdI7vlr7skrmeefzHNn+n8bl/lPC2dQOEXrEPCMtfYtY8w64NURxsxjcMWrCjg9dF8J8DmQNdxhQ2vtywxutB8XxpiD1tpl4/V6sSaZ55/Mc4fknr/mnpxzh+SefzLPHWJr/uEcR9sL5Bpj9gEHhtkID7AZeHho7CvW2gEGC9T/wWD5+q+RiSsiIiISP+66omUHNxY9f/O2MeZBYBtQYa09MzSmH9h02+O6GCxfIiIiIklpLDvDDwKL+eonD+PBuB2mjFHJPP9knjsk9/w19+SVzPNP5rlDDM3f3Osn4URERERkeGM714GIiIiI3FXSFC0z6BfGmMPGmJ85nWc8GWPmG2NajTHVQ1+Vd39U/DPGPGeMuWSMSUu2n/9tc0+an78xJtUY8+uheb4zdDuZfu63z39REv3sM40xHwzNc5sxptgY8zNjzCFjzC+HLieXsIaZf4kx5sKXfvZrnc4YbcaYbxtj2mPt/T5piha3XkroBWPM8BdOSkyZwL9Ya9cMfZ12OtA46QVunpst2X7+X557Mv38vwP/f3t3DJJVFIZx/P8QbQUhRDWFS00h6GCRqEtEIeQUBS0RQVGTNDgGQVtQQzQ2tAQSSQhBg5UR2SAptbS0JNQQ7ZH1NpzzkVw/ouV+9+Oe5zd5tvPy+B1fj9f7shoRY6Q39s5SVu7V+qcpJ/tdwGyufRE4DkxExAgwBIw2ubkeqNY/BjzdlP2rZrdXL0l7gbPAZ/rsvC+p0eo2SqgUO4ATkl5Lms9DwFsvIhaAn3lZVP6V2kvKf400YxVgAwgKyp2t9W9QSPYRsQ78kLQCHAO+Ay8l7SH90jHS5P7q1qX+VWBY0pKkRUkDze6wdreAGeA3fXbel9RodUYC/WuUUFt9I73f7CjpID7d8H6a4PwLyD8i3kfEuqRpUsa/KCj3LvUvUEj2ABHxMd9gPQcGSNkPAQ9pefawpf5RYC4ixoEHwJVGN1cjSReBZ51XTtFn531JjVZnlNAnYCddRgK1VUSsRcS9vHxLGo1UGueftD5/SedIfza5TIG5b66/pOwlDUranZdPgOuk242vpFvdVmffpf4jEXEjr1udPXAKOC/pBWlKzQx99LkvqdFaAUbzMwrjwLuG99Mz+cHoS3l5mL/P7pTE+U/mYEMAAADaSURBVCetzl/SftLIsGv5ZctF5V6tv6TsSTc4V/PXw8B9UsP5AZggfS+0WbX+L5Ju5nWrs4+IqYiYjIhJ0lzmKfroc19So/U/o4Ta6jFwUtIycBCYa3g/TXD+ZeR/ATjU+U8r4ABl5V6tfzvlZP8IGJS0RPpBe5s0a3cZ2Ea61Wmzav13gH2S3gBngLtNbq7H+uq89wtLzczMzGpS0o2WmZmZWU+50TIzMzOriRstMzMzs5q40TIzMzOriRstMzMzs5q40TIzMzOriRstMzMzs5q40TIzMzOryR9hy9UgyaH1IAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "xs = np.linspace(0, 40, 100)\n",
    "rv = stats.chi2(df=n-1)\n",
    "ax.plot(xs, rv.pdf(xs), color='gray')\n",
    "hist, _, _ = ax.hist(sample_y, bins=100,\n",
    "                     range=(0, 40), density=True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:49.886029Z",
     "start_time": "2018-08-14T18:04:49.881325Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(91.525, 337.596)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = stats.chi2(df=n-1)\n",
    "lcl = (n-1) * u_var / rv.isf(0.025)\n",
    "hcl = (n-1) * u_var / rv.isf(0.975)\n",
    "\n",
    "lcl, hcl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:50.027045Z",
     "start_time": "2018-08-14T18:04:49.887014Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAJCCAYAAAAP/PnVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XtspPt91/HPr95JjpLCJKItl2PIQiMuKlNSd4TKrVg8/ME1AgpIaZBwFTiIi+D8UQESAtFIRfyDMAapVQLICMJF5XaYSoCEaSgpbFpv5TDQUqnQbWRCEangoQQmmdg//tjdsF777O7Px+NnvH69pNXu/PbxPN9xVsl7fs8zTqm1BgCAF/dlQw8AAHDTCCgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARndWfYKv+IqvqHfv3l31aXjJfOELX8jb3va2occA4Ja5f//+Z2utX/m841YeUHfv3s3h4eGqT8NL5sGDBxHeAFy3UsqPv8hxLuEBADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACN7gw9ADxte3s7i8Ui9+7dG3oUALiQHSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBo9NyAKqW8vZTyt0spnyil/JNHjz9aSvnBUsqHr2NIAIB1cucFjvldSY5qrd9cSvkLSf5kkleSTJN8ppTynbXWz6xyyNtmPp/n4OAgfd9nPB6n67pMJpOhxwIAHnmRgPpUko8/+vMXk9Qk35Pka5N8IsnXJRFQV2Q+n2c2m2W5XCZJ+r7PbDZLEhEFAGviuQFVa50nSSnltyUZJ/lfST6b5KvyMKzGK5xv5fb394ce4Yzj4+OcnJycWVsul3njjTdy//79gaa6Xu973/tyenq6dv/ZXGRnZ2foEQAYwIvsQKWU8nuSvC/JH0zyh5O8PckPJfn6JD92wfGvJXktSV599dU8ePDgisa9eovFYugRzng6np5cX7dZV+X09DS11hvxetf53zYAq1Nqrc8+oJT3JNmttf72R49/RZJvqrV+aynlnyb50LPugZpOp/Xw8PAqZ36p7e7upu/7c+vj8Tivv/76ABNdv+3t7SwWi9y7d2/oUQC4ZUop92ut0+cd9yI/xuBDSSaPPoX3iSS/MMm7SimfTPIDbiC/Wl3XZTQanVkbjUbpum6giQCAp73IPVB/JsmfeWr5b6xmHB7fKO5TeACwvl7oHiiu12QyEUwAsMb8JHIAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGh0Z+gBWH/z+TwHBwfp+z7j8Thd12UymQw9FgAMRkDxTPP5PLPZLMvlMknS931ms1mSiCgAbq1Sa13pCabTaT08PFzpOdbN/v7+0CNcmePj45ycnJxb39jYyObm5krOeXR0lNPT02xtba3k+RnWzs7O0CMAvKlSyv1a6/R5x7kHime6KJ6etQ4At4FLeCvwMr3D3t3dTd/359bH4/HKXuf29nYWi0X29vZW8vwA8FbZgeKZuq7LaDQ6szYajdJ13UATAcDw7EDxTI9vFPcpPAD4/wQUzzWZTAQTADzBJTwAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGt0ZegBul/l8noODg/R9n/F4nK7rMplMhh4LAJoIKK7NfD7PbDbLcrlMkvR9n9lsliQiCoAbpdRaV3qC6XRaDw8PV3qOm2x/f3/oEa7N8fFxTk5Ozq1vbGxkc3PzS4+Pjo5yenqara2t6xwPvmRnZ2foEYCBlFLu11qnzzvOPVBcm4vi6VnrALCuXMIb2G16p7u7u5u+78+tj8fjM9+H7e3tLBaL7O3tXeN0APDi7EBxbbquy2g0OrM2Go3Sdd1AEwHA5diB4to8vlHcp/AAuOkEFNdqMpkIJgBuPJfwAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGh0Z+gBGMZ8Ps/BwUH6vs94PE7XdZlMJkOPBQA3goC6hebzeWazWZbLZZKk7/vMZrMkEVEA8AJKrXWlJ5hOp/Xw8HCl5xjK/v7+0CNcyvHxcU5OTs6tb2xsZHNzc4CJzjo6Osrp6Wm2traGHuXW2tnZGXoEgEGUUu7XWqfPO849ULfQRfH0rHUA4CyX8N6Cm/oufXd3N33fn1sfj8dr8Zq2t7ezWCyyt7c39CgAcCE7ULdQ13UZjUZn1kajUbquG2giALhZ7EDdQo9vFPcpPAC4HAF1S00mE8EEAJfkEh4AQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADS6M/QAJPP5PAcHB+n7PuPxOF3XZTKZDD0WAPAmBNTA5vN5ZrNZlstlkqTv+8xmsyQRUQCwpkqtdaUnmE6n9fDwcKXnaLW/vz/0CF9yfHyck5OTc+sbGxvZ3NwcYKLhHR0d5fT0NFtbW4POsbOzM+j5Abh+pZT7tdbp845zD9TALoqnZ60DAMO7lZfw1mlnYXd3N33fn1sfj8drNed12t7ezmKxyN7e3tCjAMCF7EANrOu6jEajM2uj0Shd1w00EQDwPLdyB2qdPL5R3KfwAODmEFBrYDKZCCYAuEFcwgMAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEbPDahSygdKKZ8tpbxSSvmaUsqnSymfePTrvdcxJADAOnmRHaifSvIjj/785Un+eq31Vz/69aOrG209fexjyd27yZd92cPfP/axoScCAK7bcwOq1vrdSZaPHn55kt9YSvm+Uso/LqW8faXTrZmPfSx57bXkx388qfXh76+9JqIA4La503j8Z5Ps11q/o5TybUl+d5K/efVjtdveXv057t1LPv/5s2v/5/8kH/pQ8tGPrv78N93HPz70BABwNZoCqtb6qSSfevTwk0l++UXHlVJeS/Jakrz66qt58ODBWxjxxSwWP2vl5/j859+epFywXrNYfP78F3DGgwc/8ULHLRaLLJfLa/l3AwCX0RRQpZQPJBnXWr8zyTck+eGLjqu1fiTJR5JkOp3Wu3fvvsUxn+/evZWfInfvPrxs97T3vKfk3r1XVj/AjXf3hY565ZWH38vr+HcDAJfR+mMM/lGS31RKuZfkFyX5rqsfaX19+7cn73jH2bV3vOPhOgBwe7zQDlStdfuJh+9fzSjr74MffPj7n/pTyac/nfy8n/cwnh6vAwC3Q+tN5LfeBz8omADgtvOTyAEAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoNGdoQfg5TKfz3NwcJC+7zMej9N1XSaTydBjAcCVElBcmfl8ntlsluVymSTp+z6z2SxJRBQAL5VSa13pCabTaT08PFzpOW6K/f39oUdYqePj45ycnJxb39jYyObm5gs/z9HRUU5PT7O1tXWV43EL7OzsDD0CcMOVUu7XWqfPO849UFyZi+LpWesAcFO5hHeNXvZ3x7u7u+n7/tz6eDxueu3b29tZLBbZ29u7wukA4OrYgeLKdF2X0Wh0Zm00GqXruoEmAoDVsAPFlXl8o7hP4QHwshNQXKnJZCKYAHjpuYQHANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQKM7Qw8ALebzeQ4ODtL3fcbjcbquy2QyGXosAG4ZAcWNMZ/PM5vNslwukyR932c2myWJiALgWpVa60pPMJ1O6+Hh4UrP8bLb398feoRrdXR0lNPT02xtbZ1ZPz4+zsnJybnjNzY2srm5eV3jwbXY2dkZegS4lUop92ut0+cd5x4oboyL4ulZ6wCwKi7h3QC37Z3o9vZ2FotF9vb2zqzv7u6m7/tzx4/H41v3PQJgWHaguDG6rstoNDqzNhqN0nXdQBMBcFvZgeLGeHyjuE/hATA0AcWNMplMBBMAg3MJDwCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGt0ZegC47ebzeQ4ODtL3fcbjcbquy2QyGXosAJ5BQMGA5vN5ZrNZlstlkqTv+8xmsyQRUQBrrNRaV3qC6XRaDw8PV3oOVmd/f//az3l0dJTT09NsbW1d+7mv2/HxcU5OTs6tb2xsZHNzc4CJgKu2s7Mz9Ag0KKXcr7VOn3ece6BgQBfF07PWAVgPLuHxTEO8c9re3s5iscje3t61n/u67e7upu/7c+vj8di7VoA1ZgcKBtR1XUaj0Zm10WiUrusGmgiAF2EHCgb0+EZxn8IDuFkEFAxsMpkIJoAbxiU8AIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABrdGXoAXj7z+TwHBwfp+z7j8Thd12UymQw9FgBcGQHFlZrP55nNZlkul0mSvu8zm82SREQB8NIotdaVnmA6ndbDw8OVnuMm2d/fH3qElTo+Ps7Jycm59Y2NjWxubr7QcxwdHeX09DRbW1tXPR63yM7OztAjADdQKeV+rXX6vOPcA8WVuiienrUOADeRS3jX7GV/V7y7u5u+78+tj8fjF37t29vbWSwW2dvbu+LpAOBq2IHiSnVdl9FodGZtNBql67qBJgKAq2cHiiv1+EZxn8ID4GUmoLhyk8lEMAHwUnMJDwCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGj0QgFVSvlAKeWzpZRXykMfLaX8YCnlw6seEABg3dx5weN+KsmPPPrzNyR5Jck0yWdKKd9Za/3MKobj5TCfz3NwcJC+7zMej9N1XSaTydBjAcClvdAOVK31u5MsHz38+iTfk+Rrk3wiydetZjReBvP5PLPZLH3fJ0n6vs9sNst8Ph94MgC4vBfdgXrSOMmnk3xVko8/ekyj/f39oUe4FsfHxzk5OTmztlwu88Ybb+T+/fsXfs373ve+nJ6e3prvEcPY2dkZegTgBrtMQPVJ3p7kh/JwN+rHnj6glPJakteS5NVXX82DBw/ewogvp8ViMfQI1+LpeHpy/c2+B6enp6m13prvEcPw30vAW1FqrS92YCkfT/Ib8vCS3TfVWr+1lPJPk3zoWfdATafTenh4eBWzcgPt7u5+6fLdk8bjcV5//fULv2Z7ezuLxSL37t1b9XgAcEYp5X6tdfq84y7zYwzuJXlXKeWTSX7ADeQ8S9d1GY1GZ9ZGo1G6rhtoIgB46174El6tdfuJh7/v6kfhZfT403Y+hQfAy+Qy90BBk8lkIpgAeKn4SeQAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANLoz9AC32Xw+z8HBQfq+z3g8Ttd1mUwmQ48FADyHgBrIfD7PbDbLcrlMkvR9n9lsliQiCgDWXKm1rvQE0+m0Hh4ervQcLfb394ceIUlyfHyck5OTc+sbGxvZ3NwcYKL1cXR0lNPT02xtbV3reXd2dq71fACsn1LK/Vrr9HnHuQdqIBfF07PWAYD1cesu4a3LLsPu7m76vj+3Ph6P12bGoWxvb2exWGRvb2/oUQDgQnagBtJ1XUaj0Zm10WiUrusGmggAeFG3bgdqXTy+Udyn8ADg5hFQA5pMJoIJAG4gl/AAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKDRnaEHAC5nPp/n4OAgfd9nPB6n67pMJpOhxwK4FQQU3EDz+Tyz2SzL5TJJ0vd9ZrNZkogogGtQaq0rPcF0Oq2Hh4crPQfD2d/fv/LnPDo6yunpaba2tq78uV8Wx8fHOTk5Obe+sbGRzc3NASYC1sXOzs7QI9xopZT7tdbp845zDxTcQBfF07PWAbhaLuHxlqzinc729nYWi0X29vau/LlfFru7u+n7/tz6eDz27hPgGtiBghuo67qMRqMza6PRKF3XDTQRwO1iBwpuoMc3ivsUHsAwBBTcUJPJRDABDMQlPACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAa3Rl6AGC9zefzHBwcpO/7jMfjdF2XyWQy9FgAgxJQwJuaz+eZzWZZLpdJkr7vM5vNkkREAbdaqbWu9ATT6bQeHh6u9BzcTPv7+xeuHx0d5fT0NFtbW9c7EOccHx/n5OTk3PrGxkY2NzcHmAi4zXZ2dlZ+jlLK/Vrr9HnHuQcKeFMXxdOz1gFuC5fwGMybvZPY3t7OYrHI3t7e9Q7EObu7u+n7/tz6eDy+lneCAOvKDhTwprquy2g0OrM2Go3Sdd1AEwGsBztQwJt6fKO4T+EBnCWggGeaTCaCCeApLuEBADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAECjO0MPcFvN5/McHByk7/uMx+N0XZfJZDL0WADACxBQA5jP55nNZlkul0mSvu8zm82SREQBwA1Qaq0rPcF0Oq2Hh4crPUeL/f39oUfI8fFxTk5Ozq1vbGxkc3NzgInWy9HRUU5PT7O1tXUt59vZ2bmW8wCw/kop92ut0+cd5x6oAVwUT89aBwDWy627hLcOuw27u7vp+/7c+ng8Xov5hra9vZ3FYpG9vb2hRwGAC9mBGkDXdRmNRmfWRqNRuq4baCIAoMWt24FaB49vFPcpPAC4mQTUQCaTiWACgBvKJTwAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGh0qYAqpXxNKeXTpZRPPPr13qseDABgXd255Nd9eZK/Xmv9s1c4C1yr+Xyeg4OD9H2f8XicrusymUyGHguAG+Cyl/C+PMlvLKV8XynlH5dS3n6VQ8GqzefzzGaz9H2fJOn7PrPZLPP5fODJALgJSq21/YtK+WVJfmWt9TtKKd+W5EdrrX/zomOn02k9PDx8i2Pyovb394ce4S07OjrK6elptra2VnaO4+PjnJycnFvf2NjI5ubmys4LL4udnZ2hR4CVKKXcr7VOn3fcpS7h1Vo/leRTjx5+Mskvf+rkryV5LUleffXVPHjw4DKn4RIWi8XQI7xlp6enqbWu9LVcFE+P11+G7yGsmv9e57a77A7UB5KMa63fWUr5cJIfrrX+nYuOtQNFq+3t7SwWi9y7d29l59jd3f3S5bsnjcfjvP766ys7LwDr7UV3oC57D9Q/SvKbSin3kvyiJN91yeeBQXRdl9FodGZtNBql67qBJgLgJrnsJbxFkvdf8SxwbR5/2s6n8AC4jMv+GAO48SaTiWAC4FL8JHIAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGh0Z+gB4Caaz+c5ODhI3/cZj8fpui6TyWTosQC4JgIKGs3n88xmsyyXyyRJ3/eZzWZJIqIAbolSa13pCabTaT08PFzpOWizv78/9AjPdHR0lNPT02xtbQ09yoWOj49zcnJybn1jYyObm5sDTAS3187OztAj8JIppdyvtU6fd5x7oKDRRfH0rHUAXj4u4d1C6/6ObXt7O4vFInt7e0OPcqHd3d30fX9ufTwer/33FoCrYQcKGnVdl9FodGZtNBql67qBJgLgutmBgkaPbxT3KTyA20tAwSVMJhPBBHCLuYQHANBIQAEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQKM7Qw8AF1ksFtnd3U3f9xmPx+m6LpPJZOixACCJHSjW0Oc+97n0fZ++75Mkfd9nNptlPp8PPBkAPGQHak3t7+8PPcJgptNpaq1n1pbLZd54443cv39/oKngrdnZ2Rl6BOAK2YFi7TwdT4+dnJxc8yQAcLHyZv9jdVWm02k9PDxc6Tl4ubz3ve/NF7/4xXzLt3zLmfXxeJzXX399oKkAuA1KKfdrrdPnHWcHirXz7ne/O6WUM2uj0Shd1w00EQCc5R4o1s473/nOLJfLjMdjn8IDYC0JKNbSK6+84nIdAGvLJTwAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGt0ZegCu13w+z8HBQfq+z3g8Ttd1mUwmQ48FADeKgLpF5vN5ZrNZlstlkqTv+8xmsyQRUQDQoNRaV3qC6XRaDw8PV3qO67a/vz/0CJdyfHyck5OTc+sbGxvZ3NwcYKKLHR0d5fT0NFtbW0OPcivs7OwMPQLA2iil3K+1Tp93nHugbpGL4ulZ6wDAxVzCu4Sb+o59d3c3fd+fWx+Px2v1mra3t7NYLLK3tzf0KABwITtQt0jXdRmNRmfWRqNRuq4baCIAuJnsQN0ij28U9yk8AHhrBNQtM5lMBBMAvEUu4QEANBJQAACNBBQAQCMBBQDQSEABADQSUAAAjQQUAEAjAQUA0EhAAQA0ElAAAI0EFABAIwEFANBIQAEANBJQAACNBBQAQKM7Qw8A62w+n+fg4CB932c8Hqfrukwmk6HHAmBgAgrexHw+z2w2y3K5TJL0fZ/ZbJYkIgrgliu11pWeYDqd1sPDw5Weg2fb398feoQmR0dHOT09zdbW1qBzHB8f5+Tk5Nz6xsZGNjc3B5gIXl47OztDjwBJklLK/Vrr9HnHuQcK3sRF8fSsdQBuD5fwboGb9s5ue3s7i8Uie3t7g86xu7ubvu/PrY/H4xv3PQXgatmBgjfRdV1Go9GZtdFolK7rBpoIgHVhBwrexOMbxX0KD4CnCSh4hslkIpgAOMclPACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaHSpgCoPfbSU8oOllA9f9VAAAOvssjtQ35DklSTTJK+VUn7O1Y0EALDeLhtQX5/ke5J8bZJPJPm6K5sIAGDNXTagxkk+m+Srknz80WMAgFvhziW/rk/y9iQ/lIe7UT/25F+WUl5L8lqSvPrqq3nw4MFbGJHbZn9/Pz/5kz/p3w0Aa+uyAXU/yTfVWr+rlPKNSf7Gk39Za/1Iko8kyXQ6rXfv3n1LQ3I7+XcDwLq67CW8e0neVUr5ZJIfqLV+5gpnAgBYa5fagaq11iS/74pnAQC4EfwgTQCARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGgkoAIBGAgoAoJGAAgBoJKAAABoJKACARgIKAKCRgAIAaCSgAAAaCSgAgEYCCgCgkYACAGgkoAAAGpVa62pPUMp/T/LjKz0JL6OvSPLZoYcA4NZ5T631K5930MoDCi6jlHJYa50OPQcAXMQlPACARgIKAKCRgGJdfWToAQDgzbgHCgCgkR0oAIBGAoq1UUp51+PfSyk/beh5AODNuITHypVSfiLJv39q+auTfHOt9d8+OuZOku9L8o1Jfk2S31Zr/SNPPMfbkvz6JL8iyduT/PMk/zrJX661/oGVvwgAeIIdKFbq0a7Sv0nyO5J8LskfqrX++iR/L8nyiUP/UJK/W2v9fK31XyR5Zynl/U/8/U9P8r4kPz/Je5J83aO1rpTy8Ue/fu3qXxEACChW72cn+Wyt9X8l+WNJ/mIp5ZUkoyRfTJJSytcn+S1J/vITX/dHknxrKeU3P3r81UneleSdSd6R5N1Jfl2Sv5bk9yb5d7XWf7X6lwMALuGxQqWUnSR/NA//b1n+81N//d4k/zXJP0jyO5N8Pmd3pJLkR5L84iR/MslP5OHO05N+dx7G1F9LslVr/UtXOD4AvCkBxUqVUj6Yh5FzkOS3JvkrtdaTUspfTfLtSb6QZKPW+ulSyi9N8kdrra+VUn5mHu5IfXOt9fFO1T9L8sqjp/6pJP8xyS9J8v1JvqfW+q+v87UBcHu5hMeq/dwk/6nW+p+T/Ickf6+UspGH9y/9VK31v9RaP/3o2K0k80d//ulJ/vfjeHrkc0n+/qNfn0/yx5N8dx5ewvu+lb8SAHjkztAD8NL7BUm+qZTyuUePN/JwN+prk/zDUsqfePxJvCS/P8kHHv35ZyT5n0891zvz8HJf8jC+ainle5N8OMnPSXK8otcAAGcIKFbtvUl+Va31C08ullJ+oNb6jU88/rYk31trfRxBX53kvz31XD9Ra915dPzfKqW8Iw/j6Xcm+VullPc/ulkdAFbKJTxW7VuejKdSyp8rpXx/kh99Yu2f5OHu0p8upfyTi8/LAAAAaklEQVT+UsoPJvnzST7+9JOVUt5WSvmXSf5vku9I8pdqrd+bZD/J+58+HgBWwU3kDK6U8q5a69OX6y467t211v/xxOO3Pb2zBQDXQUABADRyCQ8AoJGAAgBoJKAAABoJKACARgIKAKDR/wPG6yYnTVuaoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "rv = stats.chi2(df=n-1)\n",
    "n_samples = 20\n",
    "ax.vlines(p_var, 0, 21)\n",
    "for i in range(n_samples):\n",
    "    sample_ = samples[i]\n",
    "    u_var_ = np.var(sample_, ddof=1)\n",
    "    lcl = (n-1) * u_var_ / rv.isf(0.025)\n",
    "    ucl = (n-1) * u_var_ / rv.isf(0.975)\n",
    "    if lcl <= p_var <= ucl:\n",
    "        ax.scatter(u_var_, n_samples-i, color='gray')\n",
    "        ax.hlines(n_samples-i, lcl, ucl, 'gray')\n",
    "    else:\n",
    "        ax.scatter(u_var_, n_samples-i, color='b')\n",
    "        ax.hlines(n_samples-i, lcl, ucl, 'b')\n",
    "ax.set_xticks([p_var])\n",
    "ax.set_xticklabels(['母分散'])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.136237Z",
     "start_time": "2018-08-14T18:04:50.028148Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.964"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = stats.chi2(df=n-1)\n",
    "cnt = 0\n",
    "for sample_ in samples:\n",
    "    u_var_ = np.var(sample_, ddof=1)\n",
    "    lcl = (n-1) * u_var_ / rv.isf(0.025)\n",
    "    ucl = (n-1) * u_var_ / rv.isf(0.975)\n",
    "    if lcl <= p_var <= ucl:\n",
    "        cnt += 1\n",
    "        \n",
    "cnt / len(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 正規分布の母平均(母分散未知)の区間推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.144213Z",
     "start_time": "2018-08-14T18:04:52.137771Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(64.512, 76.288)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = stats.t(df=n-1)\n",
    "lcl = s_mean - rv.isf(0.025) * np.sqrt(u_var/n)\n",
    "ucl = s_mean - rv.isf(0.975) * np.sqrt(u_var/n)\n",
    "\n",
    "lcl, ucl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ベルヌーイ分布の母平均の区間推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.151779Z",
     "start_time": "2018-08-14T18:04:52.145503Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 1, 1, 1, 1, 1, 0, 0, 1])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "enquete_df = pd.read_csv('../data/ch10_enquete.csv')\n",
    "enquete = np.array(enquete_df['知っている'])\n",
    "n = len(enquete)\n",
    "enquete[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.155661Z",
     "start_time": "2018-08-14T18:04:52.152855Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.709"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_mean = enquete.mean()\n",
    "s_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.162410Z",
     "start_time": "2018-08-14T18:04:52.156863Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.681, 0.737)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = stats.norm()\n",
    "lcl = s_mean - rv.isf(0.025) * np.sqrt(s_mean*(1-s_mean)/n)\n",
    "ucl = s_mean - rv.isf(0.975) * np.sqrt(s_mean*(1-s_mean)/n)\n",
    "\n",
    "lcl, ucl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ポアソン分布の母平均の区間推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.168492Z",
     "start_time": "2018-08-14T18:04:52.163536Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 11,  9,  9, 18, 13,  4, 10, 10,  8])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_access_df = pd.read_csv('../data/ch10_access.csv')\n",
    "n_access = np.array(n_access_df['アクセス数'])\n",
    "n = len(n_access)\n",
    "n_access[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.173705Z",
     "start_time": "2018-08-14T18:04:52.169557Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.444"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_mean = n_access.mean()\n",
    "s_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-14T18:04:52.179685Z",
     "start_time": "2018-08-14T18:04:52.174691Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9.698, 11.191)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = stats.norm()\n",
    "lcl = s_mean - rv.isf(0.025) * np.sqrt(s_mean/n)\n",
    "ucl = s_mean - rv.isf(0.975) * np.sqrt(s_mean/n)\n",
    "\n",
    "lcl, ucl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "821px",
    "left": "0px",
    "right": "1222.58px",
    "top": "111px",
    "width": "344px"
   },
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
